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Abstract 

In  many  engineering  applications,  including  surveillance,  guidance,  or  navigation,  single  stand-alone  sen¬ 
sors  or  sensor  networks  are  used  for  collecting  information  on  time  varying  quantities  of  interest,  such  as 
kinematical  characteristics  and  measured  attributes  of  moving  or  stationary  objects  of  interest  (e.g.  maneu¬ 
vering  air  targets,  ground  moving  vehicles,  or  stationary  movers  such  as  a  rotating  antennas). 

More  strictly  speaking,  in  these  applications  the  state  vectors  of  stochastically  moving  objects  are  to  be 
estimated  from  a  series  of  sensor  data  sets,  also  called  scans  or  data  frames.  The  individual  measurements 
are  produced  by  the  sensors  at  discrete  instants  of  time,  being  referred  to  as  scan  or  frame  time,  target  revisit 
time,  or  data  innovation  time.  These  output  data  (sensor  reports,  observations,  returns,  hits,  plots)  typically 
result  from  complex  estimation  procedures  themselves  characterizing  particula  r  waveform  parameters  of  the 
received  sensor  signals  ( signal  processing). 

In  case  of  moving  point-source  objects  or  small  extended  objects,  i.e.  typical  radar  targets,  often  rela¬ 
tively  simple  statistical  models  can  be  derived  from  basic  physical  laws  describing  their  temporal  behavior 
and  thus  defining  the  underlying  dynamical  system.  In  addition,  appropriate  sensor  models  are  available 
or  can  be  constructed,  which  characterize  the  statistical  properties  of  the  produced  sensor  data  sufficiently 
correct. 

As  an  introduction  to  advanced  target  tracking  techniques  characteristic  problems  occurring  in  typical 
radar  applications  are  presented;  key  ideas  relevant  for  their  solution  are  discussed. 
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1  Discussion  of  the  Basic  Ideas 

In  many  engineering  applications,  including  surveillance,  guidance,  or  navigation,  single  stand-alone  sensors 
or  sensor  networks  are  used  for  collecting  information  on  time  varying  quantities  of  interest,  such  as  kine- 
matical  characteristics  and  measured  attributes  of  moving  or  stationary  objects  of  interest  (e.g.  maneuvering 
air  targets,  ground  moving  vehicles,  or  stationary  movers  such  as  a  rotating  antennas). 

More  strictly  speaking,  in  these  applications  the  state  vectors  of  stochastically  moving  objects  arc  to  be 
estimated  from  a  series  of  sensor  data  sets,  also  called  scans  or  data  frames.  The  individual  measurements 
are  produced  by  the  sensors  at  discrete  instants  of  time,  being  referred  to  as  scan  or  frame  time,  target  revisit 
time,  or  data  innovation  time.  These  output  data  (sensor  reports,  observations,  returns,  hits,  plots)  typically 
result  from  complex  estimation  procedures  themselves  characterizing  particular  waveform  parameters  of  the 
received  sensor  signals  (signal  processing). 

In  case  of  moving  point-source  objects  or  small  extended  objects,  i.e.  typical  radar  targets,  often  relatively 
simple  statistical  models  can  be  derived  from  basic  physical  laws  describing  their  temporal  behavior  and  thus 
defining  the  underlying  dynamical  system.  In  addition,  appropriate  sensor  models  are  available  or  can  be 
constructed,  which  characterize  the  statistical  properties  of  the  produced  sensor  data  sufficiently  correct. 

As  an  introduction  to  advanced  target  tracking  techniques  characteristic  problems  occurring  in  typical 
radar  applications  are  presented;  key  ideas  relevant  for  their  solution  are  discussed. 

1.1  Sensor  Data  Exploitation:  Tracking 

Let  us  assume  a  single  stand-alone  radar  sensor  or  a  sensor  network  (distributed  or  co-located)  producing 
measurements  which  characterize  the  kinematical  parameters  of  certain  objects  of  interest  such  as  range, 
azimuth,  or  radial  velocity  with  respect  to  the  sensors’  position.  In  addition,  certain  types  of  sensors  can  be 
considered  delivering  attribute  type  measurements,  which  provide  information  of  the  objects’  characteristic 
properties  and  thus  can  be  used  for  target  classification  or  even  identification. 

For  efficiently  exploiting  the  sensor  resources  available  as  well  as  for  gaining  information  not  directly 
given  by  the  individual  sensor  reports  themselves,  appropriate  sensor  data  exploitation  algorithms  are  re¬ 
quired.  These  techniques  for  a  “post-production  processing”  of  the  sensor  data  basically  consist  in  a  tem¬ 
poral  integration  and  a  logical  analysis  of  the  data  by  exploiting  statistical  estimation  and  data  association 
methods.  In  this  context  also  the  combination  of  the  data  with  available  background  information  (“context 
knowledge”)  is  an  important  aspect.  These  sensor  data  processing  techniques  result  in  tracks,  i.e.  estimates 
of  state  trajectories,  which  statistically  represent  the  currently  available  knowledge  of  an  object  of  interest 
along  with  its  temporal  history.  Important  parts  of  the  tracks  arc  characteristic  quality  measures,  which 
quantitatively  describe  the  reliability  or  precision  of  this  information. 

In  the  two  lectues  devoted  to  target  tracking  and  data  fusion  aspects  we  address  characteristic  target 
tracking  tasks  and  sketch  the  structure  and  potential  use  of  more  advanced  algorithms  being  relevant  for 
designing  tracking  systems.  By  selected  examples  their  potential  in  view  of  real  radar  applications  are 
demonstrated.  The  results  can  directly  be  transfered  to  stand-alone  sensors  and  measurement  fusion,  but  are 
also  the  basic  elements  for  designing  more  sophisticated  sensor  data  fusion  architectures. 

1.2  Characteristic  Problem:  Ambiguity 

In  many  practical  applications  the  tracking  problem  is  characterized  by  uncertainty  and  ambiguities,  which 
are  inherent  constituents  of  the  underlying  scenario  and  the  sensor  systems  used  for  observing  targets  in 
a  region  of  interest.  The  BAYESian  approach  being  discussed  in  the  subsequent  sections  is  a  well-suited 
methodology  for  dealing  with  those  phenomena.  More  abstractly  speaking,  BAYESian  tracking  is  essentially 
a  processing  scheme  for  dealing  with  uncertain  information  (of  a  particular  type),  which  allows  to  make 
“soft”  or  “delayed”  decisions  as  long  as  it  is  not  possible  to  form  a  unique  decision  according  to  the  particular 
data  situation  currently  to  be  dealt  with. 
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Ambiguities  can  have  many  different  causes:  The  sensors  may  produce  ambiguous  data  due  to  their  lim¬ 
ited  resolution  capabilities  or  phenomena  such  a  Doppler  blindness  in  MTI  radar  applications  (MTI:  Moving 
Target  Indicator).  Often  an  additional  source  of  ambiguities  is  the  environment  of  the  object  to  be  tracked 
itself.  There  can  well  exist  dense  object  situations,  residual  clutter  nor  being  suppressed  of  the  radar’s  clutter 
filter,  man-made  noise,  or  simply  unwanted  targets  (e.g.  birds).  A  more  indirect  type  of  ambiguities  can 
raise  from  the  properties  of  the  object  to  be  tracked,  for  example  if  it  shows  qualitatively  distinct  maneuver¬ 
ing  phases.  Finally,  the  potential  background  knowledge  to  be  used  may  imply  problem-inherent  ambiguities 
such  as  road  maps  with  their  intersections  or  tactical  rules  describing  the  over-all  behavior  of  the  objects  to 
be  tracked. 

1.3  Multiple  Objects:  Sensor  Resolution 

Due  to  the  limited  resolution  capabilities  of  every  physical  sensor,  closely-spaced  objects  moving  as  a  group 
for  a  time  will  continuously  transition  from  being  resolved  to  unresolved  and  back  again. 

As  an  example  let  us  consider  a  medium  range  radar  producing  range  and  azimuth  measurements  for  a 
target  formation  consisting  of  two  targets.  In  case  of  a  resolution  conflict  an  unresolved  radar  plot  can  be 
interpreted  as  a  measurement  of  the  group  center.  For  physical  reasons  the  resolution  in  range,  azimuth, 
and  range-rate  will  be  independent  from  each  other.  In  particular,  range  and  cross-range  resolution  differ 
significancy  in  many  radar  applications.  Therefore  the  resolution  performance  of  the  sensor  is  expected  to 
depend  strongly  on  the  current  sensor-to-group  geometry  and  the  relative  orientation  of  the  targets  within 
the  group.  The  sensor’s  resolution  capability  also  determined  by  the  particular  signal  processing  techniques 
used  and  the  random  target  fluctuations.  As  a  complete  description  is  rather  complicated,  we  have  to  look 
for  a  simplified,  but  qualitatively  correct  and  mathematically  tractable  model. 

In  any  case,  the  radar  resolution  capability  in  range  and  azimuth  is  limited  by  the  corresponding  band-  and 
beam-width.  These  radar  specific  parameters  must  explicitly  enter  into  any  processing  of  possibly  unresolved 
plots.  The  typical  size  of  resolution  cells  in  a  medium  distance  is  about  50  m  (range)  and  500  m  (cross 
range).  As  in  target  formations  the  mutual  distance  may  well  be  50  -  500  m  or  even  less,  the  limited  sensor 
resolution  is  a  real  problem  in  target  tracking  [6].  Evidently,  resolution  phenomena  will  be  observed  if  the 
range  and  angular  distances  between  the  targets  arc  small  compared  with  the  resolution  parameters.  On 
the  other  hand,  the  targets  within  the  group  are  resolvable  if  the  opposite  is  true.  Furthermore  we  expect  a 
narrow  transient  region.  A  more  quantitative  description  is  provided  by  introducing  a  probability  of  being 
unresolved  Pu  depending  on  the  sensor-to-group  geometry.  For  this  quantity  in  subsection  3.2.3  a  simple 
model  is  discussed. 

As  an  example  let  us  consider  the  simplified  situation  in  Figure  la.  A  formation  with  two  targets  is 
passing  a  radar.  We  here  consider  an  echelon  formation.  R  is  the  minimum  distance  of  the  group  center 
from  the  radar.  Figure  lb  shows  the  resulting  probability  Pu(r:  R )  parameterized  by  R  =  0,  10,  30,  60 
km  as  a  function  of  the  distance  r  between  the  formation  center  and  the  radar.  The  solid  lines  refer  to  a 
formation  approaching  the  radar  (r  <  0),  the  dashed  lines  refer  to  r  >  0.  For  R  0  both  flight  phases  differ 
substantially.  Near  R  the  probability  Pu  varies  strongly  (0.85  — >-  0. 15).  For  a  radial  flight  (R  =  0)  we  observe 
no  asymmetry  and  Pu  is  constant  over  a  wide  range  (r  »  rc ).  This  discussion  makes  evident,  that  radar 
resolution  capability  strongly  depends  on  the  underlying  sensor-to-target  geometry  and  the  relative  position 
of  the  targets  within  the  group. 

1.4  Description  by  Probability  Densities 

Many  basic  ideas  and  mathematical  techniques  relevant  to  the  design  of  tracking  systems  can  be  discussed 
in  a  unified  statistical  framework  that  essentially  makes  use  of  Bayes’  Rule.  The  general  multiple-object, 
multiple-sensor  tracking  task,  however,  is  highly  complex  and  involves  rather  sophisticated  combinatorial 
and  logical  considerations  that  are  beyond  the  scope  of  this  tutorial.  For  a  more  detailed  discussion  of  the 
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a.  scenario  b.  resolution 

Figure  1:  Resolution  (Effect  of  Sensor- to-Target  Geometry) 


problems  involved  see  [4]  and  the  literature  cited  therein.  Nevertheless,  in  many  applications  the  task  can  be 
partitioned  into  independent  sub-problems  of  (much)  less  complexity. 

In  a  BAYESian  view,  a  tracking  algorithm  is  an  iterative  updating  scheme  for  conditional  probability 
densities  p(x/\Zk)-  These  densities  represent  all  available  knowledge  on  the  kinematical  state  vectors  x/  of 
the  objects  to  be  tracked  at  discrete  instants  of  time  f/  given  both,  the  sensor  data  Zk  accumulated  up  to  some 
time  tk,  typically  the  current  scan  time,  as  well  as  all  available  a  priori  information  (sensor  characteristics, 
object  dynamics,  operating  conditions,  road  and  topographical  maps,  tactical  rules,  . . . ).  Depending  on  the 
time  ti  at  which  estimates  for  the  state  vectors  x/  arc  required,  the  related  estimation  process  is  referred  to 
as  prediction  ( t /  >  tk),  filtering  ( q  =  tk),  and  retrodiction  ( q  <  tk),  respectively  [9,  7,  17,  22].  Equation  1 
illustrates  schematically  an  iterative  process  for  calculating  the  conditional  probability  densities  p(x/\Zk)' 


prediction: 

p(xk^\Zk~l) 

dynamics  model 

p(xk\Zk~l 

road/topographical  maps 

filtering: 

p(xk\Zk~l) 

current  sensor  data 
- > 

sensor  model 

p(xk\Zk ) 

retrodiction: 

Kx/-i|Zfc) 

filtering  output 
< - - - 

dynamics  model 

p(x,\Zk). 

Under  the  conditions  previously  discussed,  the  densities  have  a  particular  formal  structure:  They  are  finite 
mixtures,  i.e.  weighted  sums  of  individual  densities  that  assume  particular  data  interpretations  and  model 
hypotheses  to  be  true.  This  structure  is  thus  a  direct  consequence  of  the  uncertain  origin  of  the  sensor 
data  and  of  the  uncertainty  related  to  the  underlying  system  dynamics.  Provided  the  densities  p(x/\Zk )  arc 
calculated  correctly,  optimal  estimators  may  be  derived  related  for  various  risk  functions  adapted  to  the 
applications. 

Evidently,  iteratively  defined  tracking  algorithms  must  be  initiated  by  appropriately  chosen  a  priori  den¬ 
sities  (track  initiation,  track  extraction  [12,  1 1]).  This  is  a  relatively  simple  task  provided  the  sensor  reports 
are  actually  valid  measurements  of  the  objects  to  be  tracked.  For  low  observable  objects,  i.e.  targets  em¬ 
bedded  in  a  high  false  return  background,  however,  more  than  a  single  frame  of  observations  arc  unusually 
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(a)  forward  iteration  (b)  backward  iteration 

Figure  2:  Scheme  of  BAYESian  Density  Iteration 


necessary  for  detecting  all  objects  of  interest  moving  in  the  sensors’  field  of  view.  By  this,  a  higher  level 
detection  process  is  defined  resulting  in  algorithms  for  multiple-frame  track  extraction  (see  subsection  3.6). 

Figure  2  provides  a  schematic  illustration  of  the  BAYESian  density  iteration  scheme.  The  probability 
densities  i \Zk~l),  p{Xk\Zk),  and  p(Xk+i \Zk+l)  resulting  from  filtering  at  the  scan  times  tk-\,  bt,  and 
tk+i,  respectively,  arc  displayed  along  with  the  predicted  density  p(Xk+2\Zk+1)  (Figure  2a,  forward  iteration). 
While  at  time  tk-i  one  sensor  report  has  been  processed,  no  report  could  be  associated  to  the  track  at  time 
tk  .  Hence  a  missing  detection  according  to  a  detection  probability  <  1  is  assumed.  As  a  consequence  of  this 
lack  of  sensor  information,  the  density  p(Xk\Zk)  is  broadened,  because  target  maneuvers  may  have  occurred. 
This  in  particular  implies  an  increased  correlation  gate  for  the  subsequent  scan  time  tk+i  ■  According  to  this 
effect,  at  time  tk+ i  three  correlating  sensor  reports  are  to  be  processed  leading  to  a  multi-modal  probability 
density.  The  multiple  modes  reflect  the  ambiguity  regarding  the  origin  of  the  sensor  data  and  characterize  also 
the  predicted  density  p(Xk+2\Zk+i ).  By  this,  the  data-driven  adaptivity  of  the  BAYESian  updating  scheme 
is  clearly  indicated.  In  Figure  2b  the  density  p(Xk+2\Zk+2)  resulting  from  processing  a  single  correlating 
report  at  t,+ 2  along  with  the  retrodicted  densities  p(Xk+\  \Zk+2),  p(Xk\Zk+2),  and  p(Xk-i\Zk+2)  are  shown. 
Evidently,  newly  available  sensor  data  significantly  improve  the  estimates  of  the  past  states. 

1.5  Discussion  of  an  Example:  Dog-fights 

Track  initiation  and  maintenance  by  processing  noise  corrupted  sensor  returns  is  by  no  means  a  trivial  task  if 
the  sensor  data  arc  of  uncertain  origin  or  if  there  exists  uncertainty  regarding  the  underlying  system  dynamics. 
With  an  example  with  real  radar  data  recorded  during  a  dog-fight  exercise  we  mainly  focus  on  four  aspects: 

1 .  Data  association  conflicts  arise  even  for  well-separated  objects  if  a  high  false  return  background  is  to 
be  taken  into  account,  which  was  not  completely  suppressed  by  clutter  filtering. 

2.  Even  in  the  absence  of  unwanted  sensor  reports,  ambiguous  correlations  between  newly  received  sen¬ 
sor  reports  and  existing  tracks  are  an  inherent  problem  for  objects  moving  closely-spaced  for  some 
time. 
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3.  Additional  problems  arise  from  sensor  returns  having  a  poor  quality,  due  to  large  measurement  errors, 
low  signal-to-noise  ratios,  or  fading  phenomena  (i.e.  successively  missing  plots),  for  instance. 

4.  Besides  that,  the  scan  rates  may  be  low  in  certain  applications,  such  as  long-range  air  surveillance. 
Furthermore,  resolution  phenomena  make  the  data  association  problem  task  even  harder. 

5.  In  a  given  mission  often  clearly  distinct  maneuvering  phases  can  be  identified,  as  even  agile  targets 
do  not  always  use  their  high  maneuvering  capability.  Nevertheless,  sudden  switches  between  the 
underlying  dynamics  models  do  occur  and  are  to  be  taken  into  account. 

Figure  3a  shows  a  radar  data  set  accumulated  over  about  50  min.  Besides  many  false  alarms  (probably  due  to 
ground  clutter)  the  data  of  two  pairs  of  interceptor  aircraft  performing  an  air  combat  exercise  were  recorded. 

The  detection  probability  is  fairly  low  (40-60%).  In  addition  rather  long  sequences  of  missed  detections 
occur  (fading  phenomena).  The  clutter  density  is  about  .002  /km2.  The  data  were  collected  from  a  rotating  S- 
band  long-range  radar.  Range  and  azimuth  information  was  used  only;  the  elevation  data  were  corrupted  and 
thus  ignored.  The  radar  is  characterized  by  the  following  parameters:  scan  period:  10  sec,  range  accuracy: 
350  ft,  bearing  accuracy:  .22°,  range  resolution:  1600  ft,  bearing  resolution:  2.4°. 

Information  on  the  real  target  position  is  crucial  for  evaluating  tracking  filters.  This  is  particularly  true 
under  conditions  where  even  trained  human  observers  seem  unable  to  assess  the  filtering  output.  Here  a 
secondary  radar  was  used:  When  primary  and  secondary  radar  produced  identical  information  (within  a 
certain  correlation  gate),  the  primary  plots  received  an  ID  number.  The  target  ID  served  for  track  assessment 
exclusively  and  was  not  used  in  the  filtering  algorithm.  The  verified  primary  plots  arc  indicated  by  green 
and  blue  dots  in  Figure  3d,e  along  with  the  final  tracking  result  obtained  after  processing  the  raw  data 
(i.e.  multiple  hypothesis  tracking  (MHT)  and  and  subsequent  retrodiction). 

Figure  3b  shows  the  underlying  hypothesis  tree  formed  by  the  tracker  in  the  first  phase  of  dog-fight  1 . 
The  yellow  dots  indicate  hypotheses  related  to  target  1,  while  the  orange  dots  refer  to  target  2.  During  the 
tracking  process  the  number  of  targets  involved  (i.e.  two)  was  assumed  to  be  known.  Right  after  the  split-off 
maneuver  the  non-maneuvering  target  has  to  be  tracked  in  presence  of  strong  clutter  interference  and  is  thus 
likely  to  be  lost  if  mono-hypothesis  tracking  algorithm  were  used.  The  leaves  of  the  resulting  hypothesis 
tree  in  Figure  3b  represent  the  knowledge  of  the  targets’  kinematical  state  at  the  present  time.  The  impact 
of  retrodiction  on  the  available  knowledge  of  the  past  target  state  is  displayed  in  3c.  The  blue  dots  indicate 
hypotheses,  which  could  be  deleted  by  exploiting  sensor  data  which  became  available  after  the  time  when 
they  were  formed.  The  statistical  weighting  factors  after  their  creation  might  well  have  been  larger  than 
the  weight  of  the  hypothesis  produced  by  processing  the  true  target  measurement.  Red  and  white  are  the 
trajectories  finally  found  by  applying  retrodiction.  Evidently,  retrodiction  does  not  improve  the  estimates 
at  present.  We  can  conclude  that  the  ambiguities  inherent  in  the  sensor  data  can  be  removed  by  using 
retrodiction  techniques  at  the  expense  of  a  certain  time  delay  of  some  sensor  scans.  Even  a  delay  of  only  two 
frames  can  significantly  improve  the  filtering  output. 

1.6  Generic  Scheme  of  a  Tracking  System 

Figure  4  provides  a  schematic  overview  of  a  generic  tracking  system  along  with  its  relation  to  the  underlying 
sensor  system.  In  the  subsequent  sections  its  basic  elements  are  being  discussed  in  greater  detail. 

After  passing  the  detector  device,  which  essentially  serves  as  a  means  of  data  rate  reduction,  the  sen¬ 
sor  signal  processing  unit  provides  estimates  of  signal  parameters  characterizing  the  waveforms  received  by 
the  sensing  hardware  (e.g.  radar  antennas).  From  these  preprocessed  estimates  sensor  reports  are  formed, 
i.e.  measured  quantities  possibly  related  to  the  objects  of  interest,  that  are  input  information  for  the  tracking 
system.  In  the  tracking  system  itself  all  sensor  data,  which  can  be  associated  to  the  already  existing  tracks, 
are  used  for  track  maintenance  (prediction,  filtering,  retrodiction).  The  remaining  non-associated  data  arc 
processed  in  order  to  establish  new  tentative  tracks  (track  initiation,  multiple  frame  track  extraction).  The 
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(a)  accumulated  radar  data 


(b)  MHT  hypothesis  tree  (c)  retrodicted  trajectory 


(d),  (e)  both  dogfights  with  verified  primary  plots 


Figure  3:  A  Typical  Dog-fight  Scenario 
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Tracking  &  Fusion  System 


Figure  4:  Generic  Scheme  of  a  Tracking  System 


plot-to-track  association  unit  thus  plays  a  key  role  in  any  multiple  target  tracking  system.  Evidently,  a  priori 
knowledge  in  terms  of  statistical  models  of  the  sensor  performance,  object  characteristics  (including  their 
dynamical  behavior),  and  the  object  environment  is  prerequisite  to  both  track  maintenance  and  track  initi¬ 
ation.  Track  confirmation/termination,  object  classification/identification,  and  fusion  of  tracks  representing 
identical  information  is  performed  in  the  track  processing  unit.  The  generic  scheme  of  a  tracking  system 
is  completed  by  a  man-machine  interface  with  displaying  and  interaction  functions.  The  available  informa¬ 
tion  on  the  sensor,  the  objects  of  interest,  and  the  environment  can  be  specified,  updated,  or  corrected  by 
direct  human  interaction  as  well  as  by  the  track  processor  itself,  e.g.  as  a  consequence  of  a  successful  object 
classification. 

2  BAYESian  Approach  to  Target  Tracking 

Following  the  spirit  of  the  preliminary  discussion  in  the  introduction  we  briefly  summarize  along  which  lines 
how  we  shall  proceed: 

•  Basis:  In  the  course  of  time  one  or  several  sensors  produce  ‘measurements’  of  one  or  several  targets  of 
interest.  The  accumulated  sensor  data  an  example  of  a  ‘time  series’.  Each  targets  is  characterized  by  its 
current  ‘state’,  a  vector  typically  consisting  of  the  current  target  position,  its  velocity,  and  acceleration. 
The  target  state  is  expected  to  change  with  time. 

•  Objective:  Learn  as  much  as  possible  about  the  individual  target  states  at  each  time  of  interest  by 
analyzing  the  ‘time  series’  created  by  the  sensor  data. 

•  Problem:  The  sensor  information  is  inaccurate,  incomplete,  and  eventually  even  ambiguous.  More¬ 
over,  the  phenomena  determining  the  targets’  temporal  evolution  are  usually  not  well-known. 

•  Approach:  Interpret  sensor  measurements  and  target  state  vectors  as  random  variables.  Describe  by 
the  corresponding  probability  density  functions  (pdf)  what  is  known  about  these  random  variables. 
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•  Solution:  Derive  iteration  formulae  for  calculating  the  probability  density  functions  of  the  state  vari¬ 
ables  and  develope  a  mechanism  for  initiating  the  iteration.  Derive  state  estimates  from  the  densities 
along  with  appropriate  quality  measures. 


2.1  Probability  Densities:  Selected  Facts 

Let  us  at  first  collect  some  facts  from  elementary  probability  theory  to  be  used  in  the  subsequent  discussion: 

1.  Information  of  a  random  variable  x  is  gained  by  integrating  the  corresponding  probability  density 
function  p(x).  Integration  over  a  volume  V ,  i.e.  fv  dxp(x),  yields  the  probability  that  the  event  ‘x  e  V’ 
occurs.  With  this  interpretation  a  pdf  must  be  non-negative,  p(x)  >  0,  and  normalized,  }  dx  p(x)  =  1. 

2.  The  expectation  of  x  is  defined  by  the  integral  E[x]  =  J  dx  x  p(x),  i.e.  by  the  ‘centroid’  of  its  pdf. 
Another  important  expectation  is  the  ‘expected  error  of  the  expectation  of  x’,  i.e.  a  quality  measure  for 
E[x].  It  is  defined  by  the  integral  ( covariance  matrix)-. 

E[(x  -  E[x])(x  -  E[x])t]  =  Jdx  (x  -  E[x])(x  -  E[x])Tp(x).  (2) 


3.  A  conditional  probability  density  p(x|y)  of  a  random  variable  x  describes  how  available  knowledge 
about  another  random  variable  y  affects  our  knowledge  on  x.  The  conditional  pdf  is  defined  by: 


t(x  |y)  = 


t(x,  y) 
p(  y) 


(3) 


with  p(x,  y)  denoting  the  joint  pdf  of  both  random  variables  x  and  y. 


4.  By  writing  the  pdf  p(x)  of  a  random  variable  x  in  form  of  a  marginal  probability  density, 

p(x)  =  \dy  p(x,  y)  =  \dy  p(x|y)  p( y),  (4) 

we  are  able  to  bring  another  random  variable  y  into  the  play,  which  might  be  related  to  x. 


5. 


By  using  BAYES  formula  we  can  calculate  how  information  on  y  affects  our  knowledge  on  x,  provided 
the  pdfs  p(y|x)  and  p(x)  are  known.  It  is  a  direct  consequence  of  the  last  two  statements  and  is  given 
by: 


T(x|y)  = 


p(y|x)p(x) 

I  dx  p(y\x)p(x) 


(5) 


6.  Precise  knowledge  that  a  random  variable  x  is  equal  to  a  certain  value  x  fits  well  into  the  description  of 
uncertainty  by  means  of  probability  densities  if  Dirac’s  ^-distributions  p(x)  =  <5(x;x)  are  considered. 
In  this  case  we  have  for  suitable  function  g  :  x  i->-  g(x)  of  x:  E[g(x)]  =  \  dx  g(x)  <5(x;  x)  =  g(x). 


7.  An  important  special  case  is  the  Gauss ian  pdf  characterized  by  a  single  maximum  concentrated 
around  x.  Let  the  quadratic  form  q(x')  =  ^(x  -  x)TC“I  (x  -  x)  be  a  measure  for  the  distance  between 
the  random  variable  x  and  the  ‘center’  x.  By  q(x)  =  const,  ellipsoids  are  defined  centered  around  x, 
whose  volume  and  orientation  are  determined  by  a  symmetric  and  positive  definite  matrix  C.  As  a 
special  pdf  decaying  with  an  increasing  distance  of  x  from  x,  let  us  consider  p(x)  =  c~q(x)  /  \  dx  c~q{x) . 
Evidently,  p(x)  is  positive  and  correctly  normalized.  After  integration  we  obtain: 

p(x)  =  Jf(x-  x,  C)  =  det[27rCH  (6) 

with  an  expectation  vector  and  an  covariance  matrix  given  by  E[x]  =  x  and  E[(x  -  x)(x  -  x)T]  =  C, 
respectively.  By  this,  the  covariance  matrix  C  has  a  simple  and  intuitive  geometrical  interpre¬ 
tation.  By  considering  ‘C  — >•  O’,  a  representation  of  the  ^-distribution  is  defined:  <5(x;x)  ‘=’ 
lirmc  o'  A/"(x;  x,  C). 
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8.  For  GAUSSian  pdfs  there  exists  an  extremely  useful  product  formula  facilitating  many  calculations: 

jV(z;  Hx,  R)  ,V(x;  y,  P)  =  W(z;  Hy,  S)  ,V(x;  y  +  W(z  -  Hy),  P  -  WSWT)  (7) 
with:  S  =  HPHt  +  R  and  W  =  PHTS-1 

This  formula  can  be  proven  by  interpreting  W(z;  Hx,  R)  A/"(x;  y,  P)  as  a  joint  density  p( z,  x)  = 
p(z|x)p(x).  It  can  be  shown  that  p(z,  x)  is  a  GAUSSian  itself: 

p(z,x)  =  ((!)■,  (H/),  (p£t“P)).  (8) 

Using  p{ z,  x)  and  the  ‘matrix  inversion  lemma’  (e.g.  the  useful  book  [10])  calculate  the  marginal 
and  conditional  densities  p( z),  p(x|z).  Due  to  p{z\x)p{x)  =  p{x\z)p{z)  the  formula  is  obtained.  An 
equivalent  version  of  the  product  formula  is: 

A/"(z;  Hx,  R)  jV(x;  y,  P)  =  jV(z;  Hy,  S)  W(x;  Q^-V  +  HTR_1z),  Q)  (9) 

with:  Q  =  (P_1  +  HTR  “H)1  . 

9.  Let  x  be  a  GAUSSian  random  variable.  The  pdf  of  y  =  a  +  Ax  with  fixed  a  and  A  is  given  by 

A/"(x;  x,  X)  -  ‘l+Ax  ■  A/"(y;  a  +  Ax,  AXAT).  (10) 

Proof:  p{ y)  =  J  dx  p(x,  y)  =  J  dx  p(y|x)  p(x)  =  J  dx  <5(y;  x)  p(x)  as  we  have  precise  knowledge  of 
y  given  x  is  known.  For  ‘D  — ►  0’  we  can  thus  write:  p( y)  =  J  dx  W(y;  a  +  Ax,  D)  A/"(x;  x,  P)  = 
A/"(y;  a  +  Ax,  AXAT  +  D)  according  to  the  product  formula  Equation  7. 

2.2  Target  Tracking:  General  Problem 

Let  us  consider  a  time  series  of  measurement  sets  Z\  =  [z], ,  z"k }  related  to  target  states  x/  at  instants  of 
time  denoted  by  q,  /  =  l, ...  ,k:  Zk  =  {Zk,nk,  Z*._| ,  . . . ,  Zj,  }  =  {Z*,  n*.  }■  The  individual 

measurements  and  the  target  states  be  described  by  vectors  z\  and  X/,  respectively.  In  general:  dim  z|  < 
dirnx/. 

The  central  question  of  target  tracking  can  be  stated  as  follows:  What  can  be  known  about  the  target  states 
X/  at  time  instants  ... ,  tk-\,  A,  4+1, . . .,  i.e.  for  the  past,  at  present,  and  in  the  future,  by  exploiting  the 
sensor  data  collected  in  the  times  series  Zk  ? 

According  to  the  approach  previously  sketched,  the  answer  to  this  question  is  given  by  the  conditional 
probability  densities  p(x.i\Zk),  which  are  to  be  calculated  iteratively.  At  present  we  confine  the  discussion 
to  the  case  /  =  k\  i.e.  we  are  interested  in  the  target  states  at  the  current  time  4.  Firstly,  an  application  of 
Bayes’  formula  yields: 


p(xk\Zk)  =  p(xk\Zk.nk.Zk~l)  = 


p(Zk,  nk\xk,  Zk  l)  p(xk\Zk  }) 

J  dxk  p(Zk,  nk\xk,  Zk~x)  p{*k\Zk-1) 


(11) 


In  many  practical  cases  we  will  have:  p{Zk,  nk\xk,  Zk~l)  =  p{Zk,nk\xk).  This  means  that  the  measurement 
set  at  time  4  is  depending  only  on  the  target  state  at  this  time  and  not  on  previous  measurements.  According 
to  Equation  1 1,  the  pdf  p(Zk,  nk\xk)  evidently  needs  to  be  known  only  up  to  a  multiplicative  constant:  With 
a  function  &(Zk,  nk\xk)  oc  p(Zk,  nk\xk)  we  obtain  the  same  result.  Functions  proportional  to  a  conditional 
probability  density  in  this  sense  are  called  likelihood  functions .  The  quantities  p(xk\Zk~l)  and  Qi.Zk ,  nk\xk) 
in  Equation  1 1  have  intuitively  clear  meanings  as  sketched  in  the  following  subsections. 
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2.2.1  Target  State  Prediction 

The  pdf  p(xk\Zk~l)  is  a  prediction  of  the  target  state  for  the  time  tk  based  on  all  the  measurements  received 
in  the  past  up  to  and  including  time  tk-\.  We  can  write  p(xk\Zk~l)  as  a  marginal  density  in  order  to  bring 
the  target  state  x^_i  at  the  previous  time  tk-i  into  play: 

p(xk\Zk~l)  =  \dxk- 1  p(xk,xk^\Zk^1)  =  Jdxfc_!  p(xk\xk  i ,  Zk~l\  p(xk-i\Zk~l\.  (12) 

target  dynamics  idea:  iteration 


In  many  practical  cases  we  can  assume  p(xk |x>t— 1 .  Zk  1 )  =  p(xk\xk-i)  (Markov  property).  Furthermore  a 
GAUSSian  Markov  dynamics  is  defined  by  a  GAUSSian  transition  density, 

p(xk\xk-i)  =  J^(xk\  ¥k\k-ixk-i,  D*|fc_i),  (13) 


with  an  evolution  matrix  F^-i  and  a  dynamics  covariance  matrix  D/{|/(_|  defining  the  underlying  target 
dynamics  model.  For  a  target  state  xk  =  (rj,  rj,  rj )T  given  by  position,  velocity,  and  acceleration  vectors  in 
three  spatial  dimensions  the  following  simple  realization  is  useful  in  many  practical  applications  [14]: 


/  1  (tk-tk-i)i  \(tk-tk~ i)2i\ 
Ffc|fc-t  =  (  o  1  (tk-tk- 1)  1  ) 

\0  O  e-(<k-<k- o/o  1/ 


=  Z2(l  -  (800) 


(14) 


with  I  =  diag[l,  1, 1],  O  =  diag[0,  0, 0].  According  to  this  simple  model,  the  acceleration  rk  is  described  by 
an  ergodic  Markov  process  with  E[  rk  ]  =  0.  The  corresponding  autocoiTelation  function  decays  exponen¬ 
tially  and  is  given  by  E[  rk  rj  ]  =  Z2  cxp[— (?/<  -  ti)/0,  ]  I,  /  <  k.  This  expression  gives  a  clear  meaning  to  the 
modeling  parameters  Z  ( acceleration  width)  and  9  ( maneuver  correlation  time). 


2.2.2  Likelihood  Function 

The  likelihood  functions  Q(Zk ,  nk\xk)  describe,  what  can  be  learned  from  the  current  sensor  output  Zk,  nk 
about  the  current  target  state  xk.  In  the  special  case  of  well-separated  targets,  perfect  detection,  no  false 
returns  we  have:  nk  =  1,  Zk  =  {zk}.  With  an  idealized  sensor  model  describing  bias-free  measurements  by 
linear  functions  \lkxk  of  the  target  state,  which  are  corrupted  by  GAUSSian  white  noise  characterized  by  a 
measurement  error  covariance  matrix  R/c,  the  likelihood  function  is  given  by: 

t?(z/tl Xk)  oc  W( zk,  Ukxk,  Rk).  (15) 

The  possibly  time-dependent  matrix  FL  is  called  measurement  matrix  and  defines,  which  characteristic  prop¬ 
erty  of  the  target  is  currently  being  measured.  At  different  times  the  quality  of  the  sensor  measurements  itself 
may  change  as  well  as  the  accuracy  R/(  by  which  they  are  measured. 


2.2.3  Combination  of  Densities 

According  to  these  considerations  we  are  in  principle  able  to  calculate  conditional  pdf  p{xk\Zk)  iteratively, 


p(xk\Zk) 


#(Zk,nk\xk)  }  dxk-\  p(xk\xk-i)  p(xk-\\Zk  ‘) 

|  dxk  #(Zk,nk\xk)  J  dxk-\  p(xk\xk^)  p{xk^\Zk-x)‘ 


(16) 


by  combining  the  following  pieces  of  evidence: 


p(xk-\\Zk~l) 
p(xk  |x*_ 1 ) 
#(Zk,nk\xk) 


available  past  knowledge 
target  dynamics  model 
sensor  data,  sensor  model. 
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2.3  A  Realization:  Kalman  Filter 

The  well-known  KALMAN  filter  is  a  straight-forward  realization  of  the  general  tracking  scheme  previously 
sketched  in  the  case  of  well-separated  targets,  a  GAUSS-MARKOV  target  dynamics,  perfect  detection,  no 

false  returns.  Hence  Zk  is  a  time  series  of  single  measurements:  Zk  =  { z | . z/(  }•  It  will  become  clear 

below,  that  in  this  context  GAUSSian  pdfs,  p(xk\Zk)  =  J\f  (xk-,xk\k ,  Vk\k),  represent  the  available  knowledge 
at  each  time  tk.  They  are  to  be  calculated  iteratively  according  to  the  following  scheme: 

2.3.1  Track  Initiation 

At  the  beginning  of  the  iteration  the  pdf  p(xo|Z°)  =  J\f(x 0;  xo|o,  Po|o)  has  to  describe  the  initial  ig¬ 
norance.  In  many  cases  this  is  possible  by  choosing  a  ‘large’  covariance  matrix  Po|o-  More  strictly 
speaking,  we  initialize  the  interation  by  xo|o  =  (zo,o,  o)T,  where  zq  denotes  the  first  measurement,  and 
Po  =  diag[R<),  (vmax)2 1,  (qmax)2 1],  respectively.  In  Po|o  the  matrix  Ro  is  the  measurement  error  covariance 
matrix  of  the  first  measurement  zo,  while  ignorance  about  the  initial  velocity  and  acceleration  is  modeled  by 
spheres,  whose  radius  is  given  by  the  maximum  speed  vmax  and  acceleration  qmax,  respectively. 

2.3.2  Prediction  Step 


r(  \  dynamics  model  -/  \ 

Jv  (Xfc-i;  Xfc— l|fc  — 1 .  Pfc-l|t-l)  - >  Jv  (Xfc;  Xfc|fe_l,  P*|*:_l  )  (17) 

*k\k-l’  Vk\k-1 

with:  xjfc|jfc_i  =  Ffc|t_ixfc_i|*_i  (18) 

Pjt|*:-1  =  lP/t— 1|A:— lFA:|A:— 1T  +  1  (19) 

These  formulae  directly  result  from  the  elementary  probability  facts  collected  in  subsection  2. 1 : 

p(xk\Zk~l)  =  \dxk_\  p(xk,xk-i\Zk~l)  =  \dxk_\  p(xk\xk-i)  p(xk-i\Zk~l)  (20) 

=  \dxk- 1  Zf  {xk\  Ffc|fc_iXfc_i,  i)  ^(Xfc-i;  xjt_i|*;_i ,  Pfc-i|fc-i)  (21) 

dynamics  model  filtering  at  tk-i 

=  A f(Xk ;  Ffc|fc-ixt_i|t_i,  Fytiyt-iPyt-iiyt-iFfciit-i7  +  Dfc|fc_i)  Jofxfc- 1  J^(xk-i',  b,  B)  . 

=  1  (normalization) 

In  the  last  step  we  made  use  of  the  product  formula  for  GAUSSians  (Equation  7). 


2.3.3  Filtering  Step 


,  current  measurement  zk  .  „  . 

A/  (xfe,  x^|^_j,  Pat|At — i )  dx  (xk,  Xk\k,  P^|^)  (22) 


sensor  model:  H^,  Rk 

with:  xk\k  =  xk\k-i +  Wk\k-i(zk -Hkxk\k-i),  W^-i  =  P*|yt— tHA:TSA:|/t_i'“1 

k\k—i’  t  —  HfcPfc|fc— iH/c  T  Rfc. 


Pjt|Jt  =  Pfc|fc— t  -  Wfc|t_iSfc|t_iWi 


(23) 


Also  these  formulae  directly  result  from  elementary  probability  reasoning  (subsection  2.1)  and  an  application 
of  the  product  formula  (Equation  7): 

p(zk\xk)  p(xk\Zk~l) 


p{xk\Zk)  =  p{xk\zk,Zk  x)  = 


J  dxk  p(zk\xk)  p(xk\Zk~l) 

A" ( Zk ,  HkXk,  R/t)  J\f (Xk,  Xfc|fe_i,  Pat|A: —  1 ) 

J dxk  JsT(zk\  ukxk,  Rk)  j^(xk-  xk |k_i,  Pfc|jt-i) 
J  - .. - '  v _ L _ _ / 


(Bayes’  rule) 


(24) 

(25) 


likelihood  function 


"V- 

prediction  at  t k 
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2.3.4  Retrodiction 


■A/" (x/;  xp,  P/|t)  ■  lll'A-"'1^"llp"1  ^(xz+i;  x/+i|fc,  P/+i|fe)  (26) 

prediction  output 

with:  x/|t  =  x/|/ +  W/|/+i(x/+i|fc -x/+i|/)  W/|/+i  =  P/i/F^^P,^, 

P/|fc  =  P/|/ +  W/|/+i(P/+i|*  -  P/+i|/)W^/+1 

These  update  equations  are  also  called  Rauch-Tung-Striebel  formulae  and  result  from  the  following 
considerations: 


p(x/\Zk)  =  \dx,+ 1  p(xi,  x/+i \Zk)  p(xi+i \Zk)  =  \dx,+ 1  p(x/|x/+i,Z/:)Kx/+i|Z/c) 

=  }dx/+i  p(x/|x/+i,  Zfc)  J*f (x/+i ;  x/+i|fe,  P,+i|t)  (28) 

^  V  ^ 

retrodiction  at  ti+\ 


with  p(x/|x/+i,  Zk)  given  by 


p(x/|x/+1,  Zfe) 


P(x/+i[x/)p(x/|Z/) 

|  dxi  p(X/+i|X;)p(x,|Z0 


^(x/+ 1;  F/+i|/X/,  D/+i|/)  A^(x/;  x/p,  P/p) 
j  dx i  Af  (x/+ 1;  F/+i|/X/,  D/+i|/)  ^ (x/;  x/p.  P/p) 

^  V  y  ^  V  ^ 

dynamics  model  filtering  at  time  i/ 


(29) 


An  application  of  the  product  formulae  in  Equation  29,  insertion  of  the  result  into  Equation  28,  and  a  second 
use  of  the  product  formulae  yields  the  retrodiction  update  formulae  in  Equation  27. 

2.3.5  Discussion 

We  discuss  some  characteristic  properties  of  Kalman  filtering  and  Rauch-Tung-Striebel  retrodiction: 

•  The  Kalman  filter  algorithm  is  linear  in  the  sensor  data,  i.e.  the  superposition  principle  is  valid. 

•  The  conditional  pdfs  arc  fully  characterized  by  the  related  expectations  and  covariance  matrices. 

•  KALMAN  filtering  corrects  predictions  by  the  difference  between  actual  and  expected  measurements. 

•  Variable  revisit  intervals  as  well  as  time-dependent  dynamics  and  sensor  models  arc  inherently  admit¬ 
ted. 

•  The  computational  effort  is  rather  small;  matrix  inversions  involved  can  often  be  performed  analyti¬ 
cally. 

•  Qualitatively  speaking,  the  retrodicted  densities  p(xj\Zk)  are  ‘sharper’  than  p(Xk\Zk)  and  p{xi\Z1)- 

•  For  retrodiction  only  expectations  and  covariance  matrices  of  the  filtering  and  predictions  arc  used. 

•  The  sensor  data  themselves  arc  processed  in  the  filtering  step  only  and  not  needed  in  retrodiction. 

•  The  information  gain  by  retrodiction  is  driven  by  the  target  dynamics:  W/p+i  =  P/pF^^P^^. 

•  There  is  no  gain  by  retrodiction  at  present:  Tracks  are  ‘smoothed’  at  the  expense  of  some  delay. 

•  Retrodiction  has  potential  applications  for  target  classification/identification/IFF  from  track  data. 
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2.3.6  Non-linearities 

In  a  radar  system  typically  target  range,  azimuth,  range-rate,  and  eventually  target  elevation  are  measured. 
These  quantities  are  easily  described  in  a  polar  coordinate  system,  while  for  modeling  the  target  dynamics 
Cartesian  coordinates  better  suited.  Therefore,  it  seems  to  be  convenient  to  perform  the  prediction  step  of 
the  density  update  in  the  Cartesian  dynamics  system  and  the  filtering  in  sensor  coordinates  according  to  the 
following  scheme: 


dynamics  system: 

=(x  ,y  ,x  ,y) 


sensor  system: 

Xs =(r  ,(p  ,r  ,cp) 


I  zk~l) 


scan  k  -  1 


p{4  K-i) 

dynamics 


p{xdk\Zk  ‘) 

P«  iz‘) 

d 

^di-s 

p(4 12fe_1) 

'  sensor 

p(4\z“) 

scan  k 

scan  k 

(30) 


The  corresponding  coordinate  transformations  ld<~s  and  ts<_d,  however,  are  non-linear  and  given  by: 


trfwtx5] 


/•  cos  cp 
r  sin  cp 

r  cos  cp—rcp  sin  cp 
r  sin  cp+rcp  cos  (p 


■ \/x2+y 2 
arctan  y/x 

(: xx+yy)/y/x2+y 2 
(xy-y±)/(x1+y2) 


(31) 


This  non-linear  character  of  the  coordinate  transformations  in  particular  implies  that  a  GAUSSian  pdf 
p(x*_il  Zk~l )  is  no  longer  a  GAUSSian  after  its  transformation  into  the  dynamics  system  and  vice  versa. 
In  order  to  circumvent  this  problem  in  ‘extended’  Kalman  filtering  the  non-linear  transformations  arc  sim¬ 
ply  linearized  by  a  first  order  Taylor  expansion  around  the  filtering  in  the  sensor  system  or  around  the 
prediction  x^_j  in  the  dynamics  system,  respectively: 

t*— j[xjp  «  trf^s[x*|A;]  +  Td<_5[x*|k]  (x*  -x*|fe)  with:  Td^s  =  did^s\xsk]k\/dxsk^k  (32) 

ts<-d[x£;]  ~  +  T(/<_j[xyt|it_1]  (xk  -  XA;|A:_1)  =  dtd<_i[xfc|yt_1]/dxA.|fc_1.  (33) 

With  this  approximation  the  ‘GAUSSianity’  of  the  densities  is  preserved  according  to  Equation  10  describing 
the  pdf  of  an  affine  transform  of  a  GAUSSian  random  variable. 

Let  us  consider  a  more  simplified  situation,  where  only  range  and  azimuth  measurements  are  available. 
These  measurements  zk  =  (rk,  cpk)T  arc  characterized  by  a  diagonal  measurement  error  covariance  matrix 
Rp  =  diagfoy,  a^]  assuming  that  range  and  azimuth  measurements  arc  independent  from  each  other.  When 
the  measurements  zk  are  transformed  into  Cartesian  coordinates,  the  corresponding  measurement  error  co- 
variance  matrix  can  approximately  be  obtained  as  follows:  Let  us  expand  t[z k]  around  the  prediction  x/t|Ar— l : 

t[z fe]  =  (rk  cos  cpk,  sin  cpk)T  «  t[xfe|fc_i]  +  Tfc|k_i  (zk  -x^-i),  (34) 

where  the  corresponding  Jacobi  matrix  can  be  written  as  the  product  of  a  rotation  and  a  dilation: 

j.  _  dt[xfc|t_!]  _  /  cos (pk\k—i  -rk\k-i^<Pk\k-l  \  _  l  cos (Pk\k—\  —  sin cpk\k-\  \  /  1  0  \ 

lk\k-\  —  dxk\k-l  ~  \  sin^|t_i  rk\k-i  cos q>k\k-l  )  ~  y  sin<j»fc|fc_i  cos<pjfc|fc_i  )  ^  0  rk\k-i  )'  ^  ’ 

rotation  dilation  Sr 


According  to  Equation  10  the  measurement  error  covariance  in  Cartesian  coordinate  is  depending  on  time 
(i.e.  on  the  predicted  target  range  rk\k-\  and  azimuth  cpk\k-\)  as  well  as  on  the  underlying  sensor-to-target 
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geometry.  It  is  given  by: 


R 


-  Tw-.R'TV,  =  DA  R-  S  rDl  =  D,  (  ’l  )  DJ. 


(36) 


As  a  direct  consequence,  the  Cartesian  measurement  error  ellipses  typically  increase  with  increasing  range. 
In  certain  applications  it  may  be  useful  to  deal  with  different  measurement  accuracies,  depending  on  the 
tracking  task  under  consideration,  such  as  search,  acquisition,  or  high-precision  tracking  for  phased-array 
radar  [4]. 

There  exist  more  advanced  methods  for  dealing  with  non-linearities  such  as  “particle  filtering”  or  “un¬ 
scented  KALMAN  update  equations  yield  in  this  case: 


Initiation: 

Xl|l 

=  Zl 

Prediction: 

Xfc|Jfc-l 

=  x;t_i|*;_i,  Pfc|fc-i  =  Pfc-i|fc-i,  k  =  2,  3, . . . 

Filtering: 

Xfc|* 

=  Xfc_l|fc_l  +  Wjfc-l  (zk  -  Xfc_l|t_l)  =  Pjt|/t(R/(1Z/t  +  PA;l1|fc_1xA;— 1|*— 1  ) 

=  PfcikZf=iRr1*i 

Pfclfc 

=  Pjt— i|fc— i  -  Wfc_i  (Pfc— i|jfc— i  +  R^)WJ_1 

=  P*-l|*-l  -  Pfc— l|fc— l(P*  — 1|*  —  1  +  Rfc)  'Pfc-llfc-l  =  (Pfciq*;-!  +  R/t1)  1 

HXl.KrT1- 

In  the  last  step  we  made  use  of  the  matrix  inversion  lemma.  From  these  considerations  it  becomes  clear  that 
in  case  of  stationary  targets  the  Kalman  filter  is  equivalent  to  a  weighted,  recursive  arithmetic  mean  of  the 
sensor  data.  The  related  error  covariance  matrix  is  a  harmonic  mean  of  the  corresponding  measurement  error 
covariance  matrices.  We  collect  some  observations: 

•  If  all  measurement  covariances  Rt,  i  =  1, . . . ,  k  are  identical,  we  observe  the  expected  ‘square-root’ 
law: 


Pk\k  =  R/k.  (37) 

•  If  all  measurement  error  ellipses  involved  differ  significantly  in  the  geometrical  orientation  relative  to 
each  other,  a  much  larger  effect  can  be  observed. 

•  The  ‘statistical  intersection’  of  error  ellipses  is  described  by  calculating  the  harmonic  mean  of  the 
related  error  covariance  matrices: 


k 

5>rV  (38) 

j=i 

•  In  the  limiting  case  of  very  narrow  measurement  error  ellipses  the  triangulation  of  the  target  position 
from  bearings  is  obtained  (—>  multiple  sensor  data  fusion). 

•  These  considerations  are  also  valid  in  3D  and  for  more  abstract  measurements. 

Let  us  consider  an  example:  The  target  position  in  Cartesian  coordinates  be  given  by  r  = 
r[cos(|),  sin(|)]T  =  r/V 2(1, 1)T.  The  measurement  error  covariances  of  sensor  1  and  2  is  given  by  Ri  = 
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D!SD[  and  R2  =  D2SDj  with  Dj  =  D(f)  =  ( \  ),  D2  =  D(^)  =  ±  (_\  } ),  S  =  diag[crr2,  (r^)2]. 

We  hence  obtain  the  ‘fused’  measurement  error  covariance 


R-1  =  R71  +R21  =  D[S "  D!  +DjS_1D2  -  D[(S_1  +DiDTS_1D2D[)Di 


=  Di  ( 


+{r<rv) 

0 


0 

1Hraq,)~ 


)Dl 


(39) 

(40) 


That  means  R  is  a  sphere  with  radius  E  given  by  ^  =  -L  +  1  Let  us  consider  the  following  special 

cases:  ‘triangulation’  (cy  »  rav)  ->  R  =  (m^)2,  ‘large  distance’  ( r  »  oy/o^)  ->■  R  =  of. 

A  practically  important  problem  is  the  following:  If  there  are  more  than  one  target  in  the  common  field  of 
view  of  both  sensors,  not  every  intersection  of  bearing  beams  actually  corresponds  to  a  real  target  position. 
For  more  details  and  possible  solutions  of  resulting  “deghosting  problem”  see  [1,2]. 


2.3.7  Expectation  Gates 

Kalman  filtering  provides  also  the  means  for  calculating  the  statistical  properties  of  expected  measure¬ 
ments.  The  corresponding  pdf  is  itself  the  basis  for  calculating  an  expectation  gate  containing  an  expected 
measurement  with  a  given  probability  ( correlation  probability  Pc).  The  conditional  pdf  of  an  expected  mea¬ 
surement  zk  at  time  tk  given  the  accumulated  sensor  data  up  to  and  including  the  time  tk-\  can  be  calculated 
by: 

p(zk\Zk~l)  =  \dxk  p(zk,xk\Zk~l)  =  \dxk  p{zk\xk)  p{xk\Zk~l)  (41) 

=  \dxk  JsT(zk\  Hkxk,  Rk\  yV(xfe;  x*|fe_i,  P*|*-t)  (42) 

V  v  J  v  V  ^ 

likelihood:  sensor  model  prediction  for  time  tk 

=  J^(zk\  HfcXfc|*;_i,  Sfc|fc_i)  with:  Sk\k-i  =  H*Pjt|fc_iHj  +  R*  (43) 

according  to  the  product  formula  (Equation  7).  Evidently,  v*|*;— 1  =  zk  —  Hkxk\k-\  is  a  GAUSSian  random 
variable  with  zero  expectation  and  the  covariance  matrix  S/q/t-i-  Being  the  difference  between  actual  and 
expected  measurement,  it  is  called  innovation,  Sjt|fc_i  is  thus  referred  to  as  innovation  covariance.  By 

llv/ci/c-ill2  =  vJA;_1S“1fe_1vfc|fc_i  <  A(PC)  (44) 

an  ellipsoid  is  defined,  which  contains  the  target  measurement  zk  with  probability  Pc.  The  actual  size  of  the 
gate  parameter  A{  Pc)  for  a  given  value  of  Pc  can  be  taken  from  a  /2 -table  [4], 


3  Elements  of  Multiple  Hypothesis  Tracking 

As  an  example  let  us  consider  6  sensor  reports  produced  by  two  closely-spaced  targets  at  time  tk  (Figure  5). 
This  single  frame  of  observations  is  by  no  means  uniquely  interpretable.  Among  other  feasible  interpre¬ 
tation  hypotheses  the  black  dots  could  be  assumed  to  represent  real  position  measurements  of  the  targets, 
while  all  other  plots  are  false  (Figure  5  a).  The  asterisks  indicate  the  predicted  target  positions  provided  by 
the  tracking  system.  Under  the  statistical  assumptions  previously  discussed,  the  expected  target  measure¬ 
ments  are  normally  distributed  about  their  predictions  with  a  covariance  matrix  Sjt|*_i  determined  by  the 
related  state  prediction  covariance  and  the  measurement  error.  As  any  prediction  uses  assumptions  on  the 
underlying  system  dynamics,  both  the  sensor  performance  and  the  dynamics  model  enter  into  the  statistics  of 
the  expected  target  measurements.  A  natural  scalar  measure  for  the  deviation  between  the  predicted  and  an 
actually  received  measurement  is  given  by  ||v/t|/t-i||2  =  vJ|Jt_1S7|1fc_1Vfc|t_i,  also  called  Mahalanobis  norm. 
Gating  means  that  only  those  sensor  returns  arc  considered  for  track  maintenance,  whose  innovations  are 
smaller  than  a  certain  predefined  threshold  A(PC). 
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(a)  two  resolved  targets  (b)  two  unresolved  targets 

Figure  5:  Sensor  Data  of  Uncertain  Origin  with  Competing  Interpretations 


Competing  with  the  previously  discussed  data  interpretation,  however,  there  exist  many  other  feasible 
association  hypotheses;  for  instance,  the  targets  could  have  produced  a  single  unresolved  measurement  as 
indicated  in  Figure  5b,  all  other  plots  being  false  returns.  Alternatively,  one  of  both  targets  might  not  have 
been  detected  or  no  target  detection  might  have  occurred  at  all,  the  gates  containing  false  returns  only.  The 
correlation  gates  and  thus  the  ambiguity  of  the  received  sensor  data  are  the  larger  the  more  false  returns  and 
missed  detections  or  even  successively  missed  detections  must  be  taken  into  account,  if  the  measurement 
errors  involved  or  the  data  innovation  intervals  are  large,  or  if  uncertainty  regarding  the  target  dynamics 
model  or  agile  targets  exists.  As  will  become  clear  below,  the  innovation  statistics  related  to  a  particular 
interpretation  hypothesis  is  essential  to  evaluating  its  statistical  weight. 

3.1  Ad-hoc  Approaches 

For  dealing  with  sensor  data  of  uncertain  origin  several  well-established  ad-hoc  methods  exist  which  are 
implemented  in  numerous  operational  tracking  systems.  Under  benign  conditions  gating  can  be  sufficient 
for  separating  real  target  measurements  from  competing  sensor  returns.  The  resulting  plot  is  then  processed 
by  Kalman  filtering  or  one  of  its  derivatives.  In  the  previous  example  (Figure  5)  two  sensor  reports  can 
be  excluded  by  this  measure.  Evidently,  the  gate  must  be  sufficiently  large,  otherwise  the  real  plot  might  be 
excluded  from  processing.  By  Nearest  Neighbor  (NN)  filters  [4]  only  the  measurement  having  the  smallest 
innovation  is  processed  via  Kalman  filtering  if  competing  returns  exist  in  the  gates.  This  approach  fails, 
however,  if  one  of  the  interpretation  hypotheses  indicated  in  Figure  5  is  true.  (Joint)  Probabilistic  Data 
Association  (PDA,  JPDA)  filters  [3]  are  adaptive  mono-hypothesis  trackers  that  show  data-driven  adaptivity 
in  case  of  data  association  conflicts  as  will  be  discussed  below. 

A  more  rigorous  Bayesian  approach,  capable  of  handling  challenging  conditions  as  sketched  in  the  intro¬ 
duction,  leads  to  Multiple  Hypothesis  Tracking  (MHT)  discussed  below  [16].  The  ad-hoc  methods  mentioned 
(KF,  NN,  PDAF,  JPDAF)  quite  naturally  prove  to  be  limiting  cases  of  this  more  general  approach. 

3.2  Sensor  Modeling 

A  statistical  description  of  what  kind  of  information  is  provided  by  the  sensor  systems  is  prerequisite  to 
processing  of  the  n^  sensor  output  data  Z/c  =  { zk  }  l  consecutively  received  at  discrete  instants  of  time  t^. 
For  the  sake  of  simplicity,  our  discussion  and  terminology  is  confined  to  point-source  objects,  small  extended 
objects,  possibly  unresolved  closely-spaced  objects,  or  small  clusters  of  such  objects;  i.e.  we  consider  “small 
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targets”  following  Oliver  Drummond’s  definition  [8].  The  underlying  statistical  models  essentially  determine 
the  feasible  interpretations  of  the  received  sensor  reports. 

3.2.1  Detection  Process 

The  detection  process  is  essentially  a  means  of  rata  rate  reduction  and  therefore  prior  to  any  further  sensor 
signal  processing  which  results  in  sensor  reports  possibly  related  to  objects  of  interest. 

For  resolved  objects  and  a  given  association  hypothesis  each  sensor  return  can  be  associated  to  exactly 
one  individual  object,  whereas  for  unresolved  closely-spaced  targets  a  given  report  may  correspond  to  several 
objects.  It  is,  thus,  reasonable  to  introduce  different  detection  probabilities  for  resolved  and  unresolved 
objects:  Pd,  Pjj .  Moreover,  the  detection  process  and  the  production  of  measurements  (fine  localization  by 
monopulse  processing,  for  instance  [5,  p.  119  ff.])  arc  assumed  to  be  statistically  independent. 

Let  FoV  (Field  of  View)  denote  a  region  large  enough  to  contain  all  relevant  sensor  reports  at  scan  k. 
False  detections  or  detections  produced  by  unwanted  objects  are  assumed  to  be  equally  distributed  in  FoV 
and  independent  from  scan  to  scan.  Moreover  let  the  number  of  false  returns  in  FoV  be  Poisson-distributed 
according  to  the  distribution 

PF{nk)  =  (|FoV|  pF )"*  e-|FoV|^,  (45) 

nk\ 

where  |FoV|  denotes  the  volume  of  FoV  and  pF  the  spatial  false  return  density. 

In  a  typical  radar  application  a  simple  quadrature  detector  decides  on  target  detection  if  the  received 
signal  strength  exceeds  a  certain  threshold:  a2k  >  Ad.  For  a  given  fluctuation  model  of  the  radar  cross 
section  of  the  targets,  the  detection  probability  depends  on  the  signal-to-noise  ratio  of  the  sensor  and  the 
detector  threshold,  while  the  false  alarm  probability  is  a  function  of  the  detector  threshold  alone.  Provided  the 
received  signal  strength  is  accessible  for  the  tracking  system,  it  can  be  used  as  input  information  for  adaptive 
threshold  control  [4],  for  discriminating  of  false  returns  [12,  13]  or  for  phased-array  energy  management 
[19].  A  more  detailed  discussion  will  be  given  in  the  second  talk  (subsection  1.2.4). 

3.2.2  Measurements 

For  a  resolved  object  let  zk  be  a  bias-free  measurement  of  its  kinematical  state  vector  xk  at  time  tk  with  an 
additive,  normally  distributed  measurement  error: 

z k  =  H kxk  +  wk,  wfe  ~  N(0,Rk )  ^  p(z/t|xfc)  =  A/"( zk\  Ukxk ,  R/c)  (46) 

as  before  in  the  case  of  KALMAN  filtering.  In  case  of  a  non-linear  relationship  between  the  target  state  and 
the  measurement,  zk  =  h^(x^),  we  thus  have  to  deal  with  data-dependent  measurement  matrices. 

In  case  of  a  resolution  conflict  we  interpret  an  unresolved  plot  zsk  at  time  tk  as  a  measurement  of  the 
group  center  [20],  i.e. 


z\  =  HgXfc  +  u  l  with  HgX£  =  ±H(x[  +  x2k),  (47) 

where  ujf  N(0,  Rg)  denotes  the  measurement  error  characterized  by  a  corresponding  group  measurement 
error  covariance  matrix  Rg  [20].  Let  H  be  the  underlying  measurement  matrix  defined  by  HxJ,  =  (r'k,  cp'k). 

Monopulse  angle  estimation  techniques  provide  approximately  bias-free,  normally  distributed  angular 
measurements.  While  in  principle  high  precision  monopulse  measurements  are  available  also  in  range,  in 
many  radar  systems  the  range  measurement  errors  are  a  supeiposition  of  errors  uniformly  distributed  in  the 
related  range  cells.  For  convenience  and  without  significant  degradation  of  the  tracking  process,  however, 
in  many  cases  range  measurement  errors  can  be  assumed  normally  distributed  with  a  standard  deviation  ar 
which  must  not  be  chosen  too  optimistically.  For  a  more  detailed  discussion  see  [4]  [chapter  2]. 


RTO-EN-SET-086 


2  - 19 


Advanced  Target  Tracking  Techniques 


ORGANIZATION 


3.2.3  Sensor  Resolution 


We  expect  that  the  resolution  performance  of  the  sensor  strongly  depends  on  the  current  sensor-to-group 
geometry  and  the  relative  orientation  of  the  targets  within  the  group.  For  physical  reasons  the  resolution  in 
range  and  azimuth  will  be  independent  from  each  other.  The  sensor’s  resolution  capability  also  depends  on 
the  particular  signal  processing  used  and  on  the  random  target  fluctuations.  As  a  complete  description  is 
rather  complicated,  we  are  looking  for  a  simplified,  but  qualitatively  correct  and  mathematically  tractable 
model. 

In  any  case  the  resolution  capability  in  range  and  azimuth  is  limited  by  the  band-  and  beam-width  of  the 
sensor  characterized  by  the  parameters  ar,  a v.  These  radar  specific  parameters  must  explicitly  enter  into  any 
processing  of  possibly  unresolved  plots. 

Resolution  phenomena  will  be  observed  if  the  range  and  angular  distances  are  small  compared  with  ar, 
av:  Ar/ar  <  1,  A cp/av  <  1.  The  targets  within  the  group  arc  resolvable  if  A r/ar  »  1,  A cp/av  »  1. 
Furthermore  we  expect  a  narrow  transient  region.  A  more  quantitative  description  is  provided  by  introducing 
a  resolution  probability  Pr  =  Pr{Ar,  Acp)  depending  on  the  sensor-to-group  geometry.  It  can  be  expressed 
by  a  corresponding  probability  of  being  unresol vabe  Pu.  Let  us  describe  Pu  by  a  Gaussian-type  function  of 
the  relative  range  and  angular  distances  [20] : 


Pr(Ar,  Acp)  =  1  -  Pu(Ar,  Acp ) 


with  Pu(Ar,  Acp)  =  exp 


log  2(  Ail)2 


exp 


-log202 


(48) 

(49) 


Evidently,  this  simple  model  for  describing  resolution  phenomena  reflects  the  previous,  more  qualitative 
discussion.  We  in  particular  observe  that  Pu  is  reduced  by  a  factor  of  2  if  Ar  is  increased  from  zero  to  ar. 
Due  to  the  Gaussian  character  of  its  dependency  on  the  state  vector  x^  the  probability  Pu  can  be  written  in 
terms  of  a  normal  density: 


Pu  =  exp  [-  log  2  [{r\  -  r2k)/ar]2\  exp  [-  log  2  [( cp\  -  cp^/a^]2]  (50) 

=  exp  j  —  log  2  (Hx^  -  Hxit2)TA-1(Hx'  -  Hx*2)]  (51) 

=  exp  [-  log  2  (Huxfc)TA_1H„xfc]  .  (52) 

Here  the  resolution  matrix  A  is  defined  by  A  =  diag(a2,  a2),  while  the  quantity  Hux/<  =  H(x[  -  x2)  can  be 

interpreted  a  measurement  matrix  for  distance  measurements. 

Up  to  a  constant  factor  the  resolution  probability  probability  P„(Xfc)  might  formally  be  interpreted  as 
the  ficticious  likelihood  function  of  a  measurement  0  of  the  distance  H(x;  -  x2k)  between  the  targets  with  a 
corresponding  ficticious  measurement  error  covariance  matrix  R„  defined  by  the  resolution  parameters  ar, 

0C(p. 

Pu(xk)  =  12/rRJ-1/2  (O;  Uuxk,  R„)  (53) 

With  R"  =  2k)g~2  =  2TfediagK2'°^l-  (54) 

According  to  a  first  order  Taylor  expansion  around  the  predicted  range  rsk^k_ l  and  azimuth  cpsk\k_l  of  the 
group  center,  the  resolution  matrix  Ac  describing  the  resolution  cells  in  Cartesian  coordinates  proves  to  be 
time  dependent  and  results  from  the  matrix  A  by  applying  a  rotation  (  around  (p"k\k_\  and  a  dilatation 

diag  [l.r*^]: 


Ac  =  D«v 


,( 


o  (r; 


k\k 


0 


)D 


T 

^lt-1  ’ 


(55) 


In  the  same  way  as  the  Cartesian  measurement  error  ellipses,  the  Cartesian  “resolution  ellipses”  depend  on 
the  target  range.  Suppose  we  have  ar  =  100  m  and  a(p  =  1°,  then  we  expect  that  the  resolution  in  a  distance 
of  50  km  is  about  100  m  (range)  and  900  m  (cross  range).  As  for  military  targets  in  a  formation  their  mutual 
distance  may  well  be  200  -  500  m  or  even  less,  resolution  is  a  real  target  tracking  problem[6]. 
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3.3  Likelihood  Functions 

The  likelihood  functions  proportional  to  the  conditional  probability  density  p(Zk,  nk\xk)  statistically  describe 
what  a  single  frame  of  n k  observations  Zk  =  {zj,  }"k_x  can  say  about  the  single/joint  state  x/<  of  the  objects 
to  be  tracked.  Due  to  the  Total  Probability  Theorem,  p(Zk,  nk  \xk)  can  be  written  as  a  sum  over  all  possible 
data  interpretations  Ek,  i.e.  over  all  hypotheses  regarding  the  origin  of  the  data  set  Zk: 

P(Zk ,  nk |x^ )  —  p(Zk,nk,Ej\xk) 
j 

=  ^  p(Zk,nk\Ej,xk)  p(Ej\xk) 
j 

As  shown  below,  the  probability  P(Ej\xk)  of  Ej  being  correct  as  well  as  the  individual  likelihood  functions 
p{Zk,nk\Ej,xk)  =  p{Zk\Ej,nk,xk )  p(nk\Ej)  directly  result  from  the  statistical  sensor  model  previously 
discussed  (eqs.  45,  46,  47,  53).  These  considerations  make  evident  that  the  determination  of  mutually 
exclusive  and  exhaustive  data  interpretations  is  prerequisite  to  sensor  data  processing.  Though  this  is  in 
general  by  no  means  a  trivial  task,  in  many  practical  cases  a  given  multiple-object  tracking  problem  can 
be  decomposed  into  independent  sub-problems  of  reduced  complexity.  We  consider  two  examples  that  arc 
practically  important,  but  can  still  be  handled  more  or  less  rigorously. 


(56) 

(57) 


3.3.1  Well-separated  Targets 

For  well-separated  objects  in  a  cluttered  environment  essentially  two  classes  of  data  interpretations  can  be 
identified  [3]: 


1.  Eq\  The  object  considered  was  not  detected,  all  nk  sensor  returns  in  Zk  are  false,  i.e.  assumed  to  be 
equally  distributed  in  FoV  (one  interpretation). 

2.  Ej,  j  =  1 , ,nk:  The  object  was  detected,  zJk  e  Zk  is  the  corresponding  measurement,  all  other 
sensor  returns  are  false  {nk  interpretation  hypotheses). 


Standard  probability  reasoning  yields: 

p(Ej\xk)  = 


1  -Pd  7  =  0 


_  f  PF(nk )  |FoV|-"‘  j  =  0 

\pF{nk  -  1)  IFoVr-1  W(z{;  Uxk,  R)  j  ±  0. 


P(Zk,  nk\Ej,  xk) 

With  ppirik )  given  by  Equation  45,  the  conditional  pdf  p(Zk,  nk |x^)  is  proportional  to  the  sum: 


nk 


p(Zk,  nk\xk)  oc  (1  -  PD)  pF  +  Pd^  7V(z{;  Hkxk,  Rk). 

7=1 


(58) 

(59) 


(60) 


up  to  a  factor  p"p  1  |FoV|  "k  e  lFoVFr  being  independent  of  the  kinematical  target  state  xk. 


3.3.2  Target  Formations 

For  a  cluster  of  two  closely-spaced  objects  moving  in  a  cluttered  environment  five  different  classes  of  data 
interpretations  exist  (xk  =  (x| ,  x~k)T)  [20]: 
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1.  Eu,  i  =  I . nk:  Both  objects  were  not  resolved  but  detected  as  a  group,  z'k  e  Zk  represents  the 

group  measurement,  all  remaining  returns  are  false  (nk  data  interpretations): 

p(Zk,nk\Eu,xk)  =  IFoVI1^  Jf (z^;  H*xfc,  R8k)  pF(nk  -  1)  (61) 

P(E„\xk)  =  £  Pu(xk)  P".  (62) 

With  P„  as  represented  in  Equation  53,  p(Zk,  nk ,  E'k\xk)  is  up  to  a  constant  factor  given  by: 

p(Zk,nt,  £>*)  oc  jV(  (<  )  ;  (  g  )*.  (  “  )  ),  (63) 

Hence  under  the  hypothesis  E‘k  two  measurements  are  to  be  processed:  the  (real)  plot  z'k  of  the  group 
center  Ugkxk  =  lH(x>  +  x2k)  and  a  (ficticious)  measurement  ‘zero’  of  the  distance  H„x/(  =  H(x[  -  x^) 
between  the  objects.  We  can  thus  speak  of  ‘negative’  sensor  information  [15],  as  the  lack  of  a  second 
target  measurement  conveys  information  on  the  target  position.  For  in  case  of  a  resolution  conflict  the 
relative  target  distance  must  be  smaller  than  the  resolution. 

2.  E0:  Both  objects  were  neither  resolved  nor  detected  as  a  group,  all  returns  in  Zk  are  thus  assumed  to 
be  false  (one  interpretation  hypothesis): 


p(Zk,nk\E0,xk )  =  Pu(xk)  (1  -  P")  pF(nk)  (64) 

PCEoN)  =  Pu(xk)  (1  -  P“).  (65) 

In  analogy  to  the  previous  considerations  we  can  write  up  to  a  constant  factor: 

p(Zk,nk,EkQ\xk)  oc  Jf  (0;  Hux,  R„).  (66) 

This  means  that  even  under  the  hypothesis  of  a  missing  unresolved  plot  at  least  a  ficticious  distance 
measurement  0  is  being  processed  with  a  measurement  error  given  by  the  sensor  resolution. 

3.  Ejj,  i,j  =  1, . . . ,  nk,  i  f  j:  Both  objects  were  resolved  and  detected,  zlk,  zJk  £  Zk  are  the  measure¬ 
ments,  nk  —  2  returns  are  false  (nk(nk  -  1)  interpretations): 

p{Zk,  nk\Ejj,  xk)  =  |FoV|2“'!<:  Jf  { z‘k\  H^x},  Rk)  (z{;  Hfcx2,  R k)  pF{nk  -2)  (67) 

P{Etj\xk)  =  (1  -  Pu(xk))  P2.  (68) 

According  to  the  factor  1  -  Pu(xk)  =  1  -  |2/rR„|“  W'(0;  Hux,  R„)  the  likelihood  function  becomes 
a  mixture,  in  which  negative  weighting  factors  can  occur.  Nevertheless  the  coefficients  sum  up  to 
one;  the  density  p{xk\Zk)  is  thus  well-defined.  This  reflects  the  fact  that  in  case  of  a  resolved  group 
the  targets  must  have  a  certain  minimum  distance  between  each  other  which  is  given  by  the  sensor 
resolution.  Otherwise  they  would  not  have  been  resolvable. 

4.  Ei0,  Em,  i  =  \ .....  nk:  Both  objects  were  resolved  but  only  one  object  was  detected,  z‘k  e  Zk  is  the 
measurement,  nk  —  1  returns  in  Zk  are  false  (2 nk  interpretations): 

p(Zk\EiQ.xk)  =  IFoVI1-"*  Jf  (z‘k,  Hfcx2,  R  k)  pF(nk  -  1)  (69) 

P(Ei0\xk)  =  E  (i  _  Pu(Xk))  Pd  (i  _  PD) .  (70) 


5.  Poo:  The  objects  were  resolved  but  not  detected,  all  nk  plots  in  Zk  are  false  (one  interpretation): 


(71) 

(72) 
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As  there  exist  +  l)2  +  1  interpretation  hypotheses,  the  ambiguity  for  even  small  clusters  of  closely- 
spaced  objects  is  much  higher  than  in  the  case  of  well-separated  objects  (nk  +  1  each).  We  thus  expect  that 
only  small  groups  can  be  handled  more  or  less  rigorously.  For  larger  clusters  (raids  of  military  aircraft,  for 
instance)  a  collective  treatment  [?]  seems  to  be  reasonable  until  the  group  splits  off  into  smaller  sub-clusters 
or  individual  objects.  Up  to  a  factor  pn£  ~  j  FoV\~nk  e~lFoV^i7  independent  of  x*  (eq.  45),  the  likelihood 
function  of  the  sensor  data. 


H/c 

p{Zk,nk\xk)  =  p(Zk,nk,E0)  +  ^  P(Zk,  Eu,  nk\xk),  (73) 

i.j= o 


is  proportional  to  a  sum  of  GAUSSians  and  a  constant: 


nk 


p(Zk,nk |x*)  oc  pl(  1  -  Pd)2(  1  -  Pu(xk))+pF(  1  -  Pfi)Pu(xk )  +  PuDpFPu(xk)  £  H [xk,  R*) 

i=  1 
nk 

+  pfPd( i  -Pd){ i  -  pu(*k))  ^  (zfc ;  R-t)+  H*x*’ R*)} 


i=  1 


+  Pjj(l  -  p„(x*))  p* (x*)  A/" (z^;  Hfcx[,  R/()  H^x2,  R*.-).  (74) 


i.j= 1 
«¥; 


3.4  MHT  Update  Equations 

The  tracking  problems  considered  here  are  inherently  ambiguous  due  to  sensor  data  of  uncertain  origin.  For 
the  sake  of  simplicity  we  concentrate  on  the  case  of  well-separated  target.  The  formalism  discussed  below, 
however,  can  directly  be  applied  to  small  target  groups  if  the  likelihood  function  in  Equation  74  is  used 
instead. 

As  in  the  examples  previously  discussed,  let  P/  denote  a  specific  interpretation  of  the  sensor  data  Z/  at 
scan  time  f/  taken  out  of  a  set  of  mutually  exclusive  and  exhaustive  interpretation  hypotheses.  Accordingly, 
the  k -tuple  P[k  =  ( Ek, ....  Pi),  consisting  of  consecutive  data  interpretations  P/,  1  <  /  <  k,  up  to  the 
time  tk,  is  a  particular  interpretation  hypothesis  regarding  the  origin  of  the  accumulated  sensor  data  Zk  = 
{Zk,  nk,  Z\ t_i ,  nk-\, . . . ,  Z \ ,  n\  \ .  Hk  is  thus  called  an  interpretation  history.  For  each  Hk  the  related  pre- 
histories  Hk-n  =  ( F/(-rt, . . . ,  Pi)  provide  possible  interpretations  of  sensor  data  Zk~n  accumulated  up  to  scan 
k  —  n.  With  llk  =  ( Ek , . . . ,  P*_n+i),  the  recent  history,  any  Hk  can  be  decomposed  in  P[k  =  (Hk,  ///<_„)■ 
Obviously,  the  each  density  p{xk\Zk)  can  be  written  as  a  sum  over  all  possible  interpretation  histories: 

P(X*|Z*)  =  ^p(xk,Hk \Zk)  =  ^ p(xk\Hk,Zk)p(Hk\Zk )•  (75) 

Hk  Hk 


p{xk\Zk)  is  thus  a  finite  mixture  density,  i.e.  a  weighted  sum  of  component  densities  p(xk\Hk,  Zk)  that 
assume  a  particular  interpretation  history  Hk  to  be  true  (given  the  data  Zk).  The  corresponding  mixing 
weights  p(Hk\Zk)  sum  up  to  one. 

3.4.1  MHT  Prediction 

Let  the  pdf  p(xk- 1 1 Zk~l)  at  time  tk-\  be  given  by  the  following  weighted  sum  of  GAUSSians: 

p(xk-i\Zk~l)  =  ^  PHk-i  ^(xt-px^,,  Ptik.fi-  (76) 

Hk -i 
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According  to  the  GAUSS-MARKOV-Model  of  the  target  dynamics  previously  introduced  (subsection  2.2.1), 
the  predicted  pdf  p(xk\Zk~l)  of  the  target  state  at  time  tk  is  given  by: 


p(*k\Zk~l) 


dxk-\  p{xk\xk~x)  p(xk~\\Zk  !)  (Markov  model) 


(77) 


=  ^  PHk_ i  Fx^^,  FPfft_1FT  +  D).  (78) 

Hk- 1 


3.4.2  MHT  Filtering 

Very  similar  to  the  proceeding  in  the  case  of  Kalman  filtering  we  obtain  the  filtering  update  equations  by 
exploiting  the  product  formula  (Equation  7): 


p(*k\Zk) 


p(Zk,nk\xk)  p{xAZk  !) 

.  -  (Bayes’  rule) 

\dxk  p(Zk,  nk\xk)  p(xk\Zk  ’) 

((1  -  Pd)  Pf  +  Pd  Tjti  ^ (< 5  Hx,,  R))  p(xk\Zk~l) 

|  dxk  ((1  -  PD)  PF  +  PD  Z-t,  (<;  Hxyt,  R))  p(xk\Zk~l) 

^  PHk  J^(Xk',  XHk,  P//J. 

Hk 


(79) 

(80) 
(81) 


The  expectation  xjjk  and  covariance  matrix  Y//k  result  from  the  KALMAN  filtering  formulae  (Equation  22), 
while  the  weighting  factors  arc  essentially  characterized  by  the  corresponding  innovations  vuk ,  the  innova¬ 
tion  covariance  matrix  Snk,  and  the  statistical  weight  of  the  corresponding  pre-history  Pnk_i : 


PHk 


PHk 


^Hk  P*Hk 


with 


P*Hk  =  PHk- 1  X 


f  (1  -  Pd)  Pf 
\Pd  JZ(,VHk\  o,  Snk) 


for 

for 


Ek 

Ek 


(82) 


3.4.3  MHT  Retrodiction 

Retrodiction  is  an  iteration  scheme  for  calculating  the  probability  densities  p(x/\Zk),  /  <  k,  that  describe  the 
past  states  x/  given  all  available  sensor  information  Zk  accumulated  up  to  a  later  scan  time  tk  >  ti,  typically 
the  current  time.  The  iteration  is  initiated  by  the  filtering  result  p{xk\Zk)  at  time  tk  and  describes  the  impact 
of  newly  available  sensor  data  on  our  knowledge  of  the  past.  In  close  analogy  to  the  previous  reasoning,  an 
application  of  the  Total  Probability  Theorem  yields: 

p{xi\Zk)  =  2  Kx/, Hk\Zk)  =  2  p(xi\Hk,Zk)  £( HkjZ^  (83) 

^ k  Hk  no  ambiguity  filtering 


The  calculation  of  p(x/\Hk,  Zk),  i.e.  given  a  particular  interpretation  history  Hk,  is  thus  as  in  subsection 


2.3.4,  i.e.  for  Pp,  =  1,  pp  =  0: 

p(x,\Hk,Zk)  =  Jf(xp,  xHk(l\k),  YHk(l\k))  (84) 

where  the  parameters  of  the  GAUSSian  are  given  by: 

xHkd\k)  =  xHkm  +  wHkm  ( xHka+i\k)  -  xHkd+i\n)  (85) 

YHkd\k)  =  YHkm  +  wHkm  (PHko+i\ k)  -  PHko+m)  wHk(i\k)T  (86) 

gain  matrix:  W Hku\k)  =  YHk(i\i)¥j+MIYHk(i+i\i)~l .  (87) 
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A  direct  consequence  of  these  considerations  is  the  notion  of  a  retrodicted  probability  [7].  Due  to  ///{  = 
(Hk ,  Hi)  and  the  Total  Probability  Theorem,  the  probability  of  Hi  being  correct  given  the  accumulated  data 
up  to  time  tk  can  be  calculated  by  summing  up  the  weighting  factors  of  all  its  descendants  at  time  tk: 


p(H,\Zk)  =  ^Jp(Hk,H,\Zk). 

Hk„ 


(88) 


This  result  can  also  be  obtained  by  considering  the  following  approximation: 


p(xi\Hk,Zk)  =  J\f  (x/;  xHk(l\k),  PHkd\k))  *  d/ (xp,  xHk(l\l),  PHjt(/|/)),  (89) 

i.e.  if  the  RTS-step  is  omitted.  This  means  in  particular:  strong  descendants  can  make  weak  ancestors 
stronger;  weak  descendants  can  weaken  also  strong  ancestors;  if  all  descendants  are  deleted,  also  the  ances¬ 
tors  die. 


3.5  Suboptimal  Realizations 

Due  to  the  uncertain  origin  of  the  sensor  data,  naively  applied  sensor  data  processing  according  to  the  pre¬ 
vious  formalism  leads  to  memory  explosions:  The  number  of  components  in  the  mixture  densities  p(xk\Zk) 
exponentially  grow  at  each  step.  Suboptimal  approximation  techniques  are  therefore  inevitable  for  any  prac¬ 
tical  realization.  Fortunately,  the  densities  resulting  from  prediction  and  filtering  are  characterized  by  a  finite 
number  of  modes  that  may  be  large  and  fluctuating  but  does  not  explosively  grow.  This  is  the  rationale  for 
adaptive  approximation  methods  that  keep  the  number  of  mixture  components  under  control  without  disturb¬ 
ing  the  density  iteration  too  seriously.  In  other  words,  the  densities  can  often  be  approximated  by  mixtures 
with  (far)  less  components.  Provided  the  relevant  features  of  the  densities  are  preserved,  the  resulting  sub¬ 
optimal  algorithms  are  expected  to  be  close  to  optimal  BAYESian  filtering. 


3.5.1  Moment  Matching 

Moment  matching  is  an  important  approximation  method,  by  which  a  pdf  p(x)  with  expectation  E;)  [x]  =  x 
and  a  covariance  matrix  Ep[(x  -  x)(x  -  x)T]  =  P  is  approximated  by  p(x)  «  jV"(x;  x,  P).  In  the  present 
context  moment  matching  is  applied  to  mixture  densities  of  the  form  p(x)  =  pu  jZ(x:  X// .  P//),  i.e.  to 
normal  mixtures.  In  this  case  x  and  P  are  given  by: 


x  =  ^  PHXH  (90) 

H  spread  term 

/  \ 

P  =  Yj  PH  {ptf  +  (XH  -  X)(XH  -  X)T  }  (91) 

H 
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(b) 


Figure  6:  Scheme  of  Moment  Matching 


due  to  the  following  calculations: 

Ep[x]  =  Jdxxp(x)  =  'YJpH\dxxJd'(x\xH,VH)  =  ^ph^h  =■  x 

H  H 

Cp[x]  =  jdxp(x)(x  -  Ep[x])(x  -  Ep[x])T  =  'YJpH\dx{x  -  x)(x  -  x)T  J^(x\xH,VH) 

H 

=  Yj  PH$dx  {(x  —  x)(x  —  x)T  -2(x-xh)(xh  -x)t}  M’ix-XH.Pn) 

H 

according  to:  \dx  (x  -  xh)(xh  -  x)T  .V(x;  xh,  P h)  =  0 

=  'YJpn\dx  {xxT  -  2 xx^j  +  xHxJj  +  xHxh  -  2xHxT  +  xxT}  A/"(x;  xH ,  PH) 

H 

=  'Yj  Pn\dx  {(x  -  xH)(x  -  xh)t  +  ( xH  -  x)(xH  -  x)T}  A^Cx;  xH,  PH) 

H 

=  Y  ph  {Ph + (xh  -  x^xh  - x)T } =  p- 

H 

Figure  6  provides  a  schematic  illustration  of  moment  matching.  A  particular'  mixture  density  p{x)  = 
cipi(x)  +  C2P2W  is  displayed  along  with  the  related  mixture  components  cipi(x),  c\p\ (x)  (Figure  6a). 
In  Figure  6b  the  mixture  p{x)  is  compared  with  the  Gaussian  density  jV"(x;  x,  P)  with  x  =  Ep  [  x  | , 
P  =  Ep  [  (x  -  x)2].  The  bars  at  the  bottom  line  indicate  the  relative  size  of  the  mixture  coefficients  ci,  C2 
in  this  example.  Evidently,  moment  matching  can  provide  a  satisfying  approximations  to  a  mixture  as  long 
as  it  is  unimodal. 

3.5.2  Single  Hypothesis 

A  radical  solution  to  the  growing  memory  problem  is  given  by  mono-hypothesis  approximations  briefly 
sketched  below: 

•  Exclusion  of  competing  sensor  data  by  testing  if  II vJfc|Jt_1 1|  >  A:  “Gating”.  If  this  is  successful  we 
obtain  Kalman  filtering  as  a  limiting  case. 
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(+)  Gating  is  very  simple  (-)  If  A  is  too  small,  the  actual  target  measurement  may  be  excluded. 

•  Forcing  a  unique  interpretation  in  case  of  conflict.  This  means  that  the  measurement  with  the  minimum 
statistical  distance  from  the  expected  measurement  is  used  for  updating:  min,-  ||vjt|fc_1||.  The  resulting 
filter  is  called  “Nearest- Neighbor- Filter  (NN)”. 

(+)  One  resultant  hypothesis.  (-)  A  hard  decision  is  taken,  which  my  be  wrong.  (-)  NN  is  not  adaptive. 

•  In  case  of  “global  combining”  all  hypotheses  arc  merged  to  one  single  representative  hypothesis.  The 
resulting  filter  is  called  “(Joint)  Probabilistic  Data  Association  Filter  (J)PADF”. 

(+)  All  sensor  data  are  used,  (+)  PDAF  is  adaptive.  (-)  Its  applicability  is  limited. 

Due  to  its  importance,  let  us  take  a  closer  look  at  the  PDAF  filter.  It  is  formally  analog  to  the  KALMAN 
filtering.  The  basic  processing  scheme  is  given  by: 

filtering  (scan  lc-1):  p(xk^\Zk~l)  =  ^/'(x*_1;  xfc_i|fc_i, Pfc_i|fc_i)  (^  initiation) 

prediction  (scan  k):  p{xk\Zk~x)  «  ^(x*;  xk\k_x,  Pa^-i)  (as  usual) 


filtering  (scan  k):  p(*k\Zk)  «  ^  pJkJ^(xk\  x^ ,  p£|/t)  «  tf(xk;  xk\k,  Pk\k) 

7=0 


where  the  quantities  xk.k,  P^,  pk  arc  to  be  calculated  as  follows: 


V 

innovation  covariance 


innovation 


gain  matrix 


weights 


With  the  combined  innovation  vk  =  2"=o  P/<  vk  we  obtain  by  moment  matching: 


(92) 


Pfc|fc  -  |fc  +  (Afclfc  xk\k)(x'k\k  *k\k)  ) 


(93) 


=  P*,*_i  -  2  +  2  P{ Wk(K  ~  vk)(vJk  ~  vfc)TWj 


(94) 


V 

spread  of  innovations 


(95) 


RTO-EN-SET-086 


2-27 


Advanced  Target  Tracking  Techniques 


ORGANIZATION 


3.5.3  Multiple  Hypotheses 

In  case  of  a  more  severe  false  return  background  or  in  a  multiple-object  tracking  task  with  correlation  gates 
overlapping  for  a  longer  time,  Bayesian  track  maintenance  inevitably  leads  to  densities  p(xk\Zk)  that  arc 
characterized  by  several  distinct  modes.  As  this  phenomenon  is  inherent  in  the  uncertain  origin  of  the  re¬ 
ceived  data,  relevant  statistical  information  would  get  lost  if  global  combining  is  applied  to  such  cases.  The 
use  of  PDA-type  filtering  is  thus  confined  to  a  relatively  restricted  area  in  parameter  space  (defined  by  pp, 
Pd,  for  instance). 

By  local  combining  of  suitably  chosen  sub-mixtures  and  pruning  of  irrelevant  mixture  components, 
however,  memory  explosions  may  be  avoided  without  destroying  the  multi-mode  structure  of  the  densities. 
Provided  this  is  carefully  done  with  data-driven  adaptivity,  all  statistically  relevant  information  may  be  pre¬ 
served  while  keeping  the  number  of  mixture  components  under  control,  i.e.  the  number  may  be  fluctuating 
and  be  even  large  in  critical  situations  but  does  not  explosively  grow.  Evidently,  PDA-type  filtering  is  a 
limiting  case  of  such  MHT-type  techniques. 

Individual  Gating.  In  a  first  step  for  avoiding  unnecessary  computational  load,  sensor  data  irrelevant  for 
a  given  track  hypothesis  are  excluded.  Individual  gating  means  that  only  those  sensor  data  are  used  for 
continuing  a  particular  track  hypothesis  H \  whose  innovations  obey:  vTH  vjjk  <  A.  The  processing 
parameter  A  must  be  tuned  to  meet  the  requirements  of  a  particular  application.  Evidently,  the  accuracy  of 
the  prediction  (depending  on  the  system  dynamics  model  and  the  previous  track  hypothesis)  and  a  priori 
information  on  the  sensor  performance  enters  into  this  decision  criterion.  Individual  gating  is  a  simple 
measure  of  pre-selecting  the  sensor  data.  It  can  be  performed  for  each  track  hypothesis  independently  before 
any  further  data  processing  takes  place. 

Pruning  Methods.  In  order  to  identify  insignificant  track  hypotheses,  first  for  each  Hk  the  weighting 
factors  pHk  are  evaluated  by  processing  the  sensor  data  within  the  gates.  This  is  done  before  the  hypothetical 
tracks  x//t,  PHk  are  computed.  Due  to  the  normalization  involved,  the  size  of  each  weighting  factor  puk 
depends  on  all  sensor  data  in  the  gates.  In  contrast  to  individual  gating,  pruning  is  therefore  applied  after 
all  weighting  factors  are  available.  In  zero-scan  pruning  track  hypotheses  are  deleted  that  are  smaller  than 
a  certain  predefined  threshold.  By  this,  an  additional  processing  parameter  is  introduced  that  must  be  tuned 
to  meet  the  requirements  of  a  particular  application.  The  limiting  case  where  the  track  hypothesis  of  highest 
statistical  weight  is  considered  only,  is  a  slightly  more  general  formulation  of  standard  Nearest-Neighbor 
filtering  as  the  hypothesis  of  a  missing  measurement  may  be  of  highest  weight.  Delayed  or  multiple-frame 
pruning  is  closely  related  to  this  procedure.  Here  we  consider  the  retrodicted  weighting  factors  for  a  past 
time  given  all  sensor  data  up  to  the  current  scan  (see  subsection  3.4.3). 

Local  Combining.  After  filtering  a  single  distinct  mode  of  p(xk\Zk)  might  be  a  superposition  of  “similar” 
mixture  components.  It  is  thus  reasonable  to  apply  local  combining  to  the  sub-mixture  producing  that  mode. 
Among  several  realizations  successive  local  combining  is  particularly  simple.  Let  us  start  with  the  mixture 
component  of  highest  statistical  weight.  In  the  order  of  decreasing  weighting  factors  a  component  is  searched 
that  is  “similar”  to  the  previous  one.  A  very  simple  scalar  criterion  for  similarity  is  provided  by: 

d(Hl’Hk)  <  K  with:  d(Hl’  Hl">  =  (xff‘  “  xiqOT(p/r‘  +  p//2)_1(xffi  -  xH2),  (96) 

where  xuu  and  P^u  denote  the  mean  and  covariance  of  the  components.  By  this,  a  third  processing 
parameter  k  is  introduced  (besides  A  and  the  pruning  parameter)  that  must  be  tuned  to  meet  the  requirements 
of  a  particular  application.  Local  combining  results  in  an  ‘effective’  component  with  an  increased  weighting 
factor.  Then  the  next  similar  component  is  searched  in  the  order  of  decreasing  weighting  factors  and  so  on. 
Having  done  this,  we  restart  the  procedure  with  the  mixture  component  having  the  second  largest  weighting 
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factor.  Due  to  the  data-driven  adaptivity  inherent  in  this  method,  MHT-type  filtering  automatically  reduces 
to  PDA-type  processing  if  PDA-processing  provides  good  approximations  to  p(x/(\Zk)- 

Objects  moving  closely-spaced  for  some  time  may  irreversibly  loose  their  identity:  When  they  dissolve 
again,  a  unique  track-to-target  association  is  impossible.  In  particular,  this  means  that  the  component  den¬ 
sities  p(x£,  H{k\Zk)  and  p(x/(,  H^\Zk)  are  nearly  identical  if  and  Ilj.  differ  only  in  a  permutation  of  the 
targets.  It  is  thus  reasonable  to  deal  with  densities  that  are  symmetric  under  permutations  of  the  individual 
targets.  By  this,  no  statistically  relevant  information  is  lost  and  the  filter  performance  remains  unchanged, 
while  the  mean  number  of  hypotheses  involved  may  be  significantly  reduced. 

3.6  Sequential  Track  Extraction 

After  solving  the  ‘track  maintenance’  problem  by  deriving  iterative  processing  schemes  for  updating  condi¬ 
tional  probability  densities,  an  important  question  is  still  open:  By  which  means  can  the  iteration  process  be 
started?  This  is  by  no  means  a  trivial  task  in  case  of  ambiguous  sensor  data. 

The  initiation  of  the  pdf-iteration  is  based  on  ‘extracted’  target  tracks,  i.e.  on  tracks  whose  existence 
is  ‘detected’  by  a  detection  process  on  a  higher  level  of  abstraction,  which  makes  use  of  sensor  detections 
accumulated  over  time.  More  strictly  speaking,  we  have  to  find  a  candidate  for  a  target  track  in  a  time  series 
of  sensor  observations  Zk  =  { Z,  }k=r  For  the  sake  of  simplicity  we  assume  for  the  time  being:  1 .  In  the  FoV 
of  the  sensors  there  is  at  most  one  object.  2.  The  sensor  data  collected  in  one  scan  arc  measured  at  the  same 
time. 

We  have  to  decide  between  two  hypotheses: 

•  hi:  Besides  false  returns,  Zk  contains  also  actual  target  measurements. 

•  ho:  There  is  no  target  existing  in  the  FoV;  all  sensor  data  in  Zk  are  false. 

Two  decision  errors  are  involved  characterizing  the  performance  of  any  test  procedure:  1.  The  conditional 
probability  that  hypothesis  h \  is  accepted  given  h\  is  actually  true  Pi  =  Prob {accept  h\\h\)  corresponding 
to  the  detection  probability  Pd  of  a  sensor.  2.  Pq  =  Prob  (accept  h\ \ho),  corresponding  to  the  false  alarm 
probability  Pp . 


3.6.1  Likelihood-ratio  Test 


We  are  looking  for  a  test  procedure  for  deciding  between  these  two  probabilities  as  quickly  as  possible 
for  given  decision  errors  Po,  P\.  Let  us  consider  the  conditional  probability  densities  p(Zk  \ho),  p{Zk\h\) 
(likelihood  functions)  and  an  intuitively  plausible  test  function  (likelihood  ratio): 


LR  (k)  = 


p(Zk\hi) 

p(Zk\hQ)' 


(97) 


Starting  from  a  time  window  of  length  k  =  1,  the  test  function  LR(k)  is  successively  calculated  and  compared 
with  two  thresholds  A  and  B: 


for  LR (k)  <  A, 
for  LR (k)  >  P, 
for  A  <  LR (k)  <  B , 


accept  the  hypothesis  ho  (i.e.  no  object  existent  in  the  FoV) 
accept  the  hypothesis  h\  (i.e.  an  object  exists  in  the  FoV) 
expect  new  data  Z/{+|  and  repeat  the  test  with  LR(k  +  1). 


This  test  procedure  ( sequential  likelihood  ratio  test)  has  the  following,  practically  important  properties: 
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1.  For  the  detection  thresholds  A,  B  and  the  decision  errors  Po,  Pi  obey  approximately  the  relationship: 


1  -  Pi  Pi 

A  « -  and  B  «  — . 

1  -  P)  Po 


(98) 


2.  The  actual  decision  time  required,  i.e.  the  amount  of  sensor  data  required,  is  a  random  quantity. 

3.  On  average  the  sequential  likelihood  test  has  a  minimal  decision  length  for  given  errors  Po,  Pi. 

4.  The  actual  choice  of  Po  (Pi)  has  impact  on  the  mean  decision  length  assuming  h\  (ho)  is  valid. 

5.  In  practice  the  parameter  Pi  is  chosen  close  to  One  for  actually  detecting  real  object  tracks. 

6.  The  parameter  Pq  should  be  small,  as  the  tracking  system  is  not  to  be  overloaded  with  false  tracks. 


3.6.2  Iterative  Updating 

For  calculating  the  likelihood  ratio,  interpretations  histories  Hk  =  { Ek,  IIk-\ }  of  the  accumulated  data 


fc-i  i 


have  to  be  considered  as  before  in  the  case  of  track  maintenance.  With  Ek  =  E V  (Object 


z*  =  \zk,z 

not  detected),  Ek  =  EJk  (z'k  e  Zk  is  the  target  measurement)  and  with  the  histories  Ilf,  we  can  write: 

_  p(Zk\hx)  _  lHk  P(Zk,Hk\hl)  _  ZHk  p(Zk\Hk,hx)p(Hk\hx) 

LR(/C)  “  p(Zk\ho)  ~  P(Zk\ho)  ~  P(Zk\ho) 

For  the  test  procedure  an  iterative  calculation  is  requested.  Standard  probability  reasoning  yields: 

f  (1  -  PD)  pF(nk)  Ek  =  El 
p(Hk\hl)  =  p(Ek\Hk-l,hl)p(Hk_1\hl)  =  p(Hk-l\h1)  {  K  „  k  k 

1  ■£  PF(nk  -  1)  Ek  =  EJk 


nk 


p(Zk\Hk,  hi)  =  p(Zk\Hk,  Zk~\  hi)  p(Zk~l \Hk_i,hi) 


=  p(Zk~l\Hk-i,  hi) 


Ek  =  E 


f  |FoV|-"‘ 

^|FoVr^+1  Af(vHk,SHk)  Ek  =  E] 


p(Zk\h0)  =  p(Zk,nk,Zk~x\ho)  =  p(Zk\nk,  Zk~x ,  h0)  p(nk\Zk~l ,  h0)  p(Zk~l\h0) 


=  IFoVr"*  PF(nk)  p(Zk~l \h0)  with:  pF(nk)  =  e-^FoV'. 

We  consider  the  following,  convenient  notation  for  multiple  sums: 


nk\ 


n\ 


h  =  (jk,  --, h )  let  us  write  ^4  =  2  "  '  2  AA-h ■ 

jk= o  ;i=0 


(99) 

(100) 

(101) 

(102) 

(103) 

(104) 

(105) 


j  k 


From  the  formulae  above  an  application  of  the  product  formulae  (Equation  7)  results  in  the  following  simple 
update  formulae  for  the  likelihood  ratio: 


Initiation: 

Updating: 


k  =  0,  j0  =  0,  Ajo  =  1 


«i+ 1 


LR(k  +  D  =  24+1=  2  2  k 

jk+ 1=0  j  k 


] . 

jk+lSk  ik 


(106) 

(107) 


j*+i 


with: 


'Jk+lik 


1  -Pd 


for:  jk+ 1  =  0 

jk+\ik  ’  ^A+ljt)  f°r:  jk+1  7^0 


(108) 
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senor  data  innovation: 

Vjk+llk  Zjk+ 1  ^jk+\^'iM\k 

(109) 

innovation  covariance: 

® ik+iik  ~  H/fc+iPjfc+iifcHy^i  +  R jk+l 

(HO) 

state  update:  xh+]]k  =  FA+1xjt 

Xj/c  —  Xjt|lt-1  4"  ^ jk)k-ivjk.jk-l 

(111) 

covariances:  Pjt+1|t  =  Fjk+lPikF]k+i 

4"  ®/'*+l  Pj*:  Pj^l  _  Wjk}k-l  ^jkik-l  ^jkjk-l 

(112) 

Kalman  gain: 

W  —  P  1 1 T  n-1 

yyjuk-i  rik\k-\  ■ 

(113) 

3.6.3  Hand-over  to  Maintenance 

The  further  proceeding  in  sequential  track  extraction  consists  of  the  following  steps: 

•  LR(&;)  is  represented  by  an  increasing  number  of  summands,  which  are  related  to  a  particular  interpre¬ 
tation  history.  The  tuple  { }  is  called  a  sub-track. 

•  For  mitigating  the  growing-memory  problem  all  approximations  are  to  be  used  which  have  been  intro¬ 
duced  for  MHT  track  maintenance,  as  far  as  they  do  no  significantly  affect  LR(/c): 

-  Individual  gating:  Exclude  data  whose  association  to  an  existing  sub-track  is  too  improbable. 

-  Pruning:  Delete  sub-tacks  which  contribute  not  significantly  to  the  likelihood  ratio. 

-  Local  combining:  Merge  similar  sub-tracks  by  using  moment  matching  according  to 

{Aj.x/.P,-},-  ->■  {2,x,P}  with:  A  =  ^  A,  (114) 

i 

P  =  j^41[P/  +  (x,-x)(.  ,.)T}.  (115) 

i  i 

•  The  test  ends  with  a  decision  in  favor  of  one  of  the  hypotheses:  ho  (no  object)  or  hi  (object  existent). 

•  After  a  track  detection,  the  pdf  for  track  maintenance  is  initiated  by  the  sub-tracks  according  to: 

A\k 

normalization  of  the  coefficients  A-jk :  p-ik  =  — — —  (116) 

An 

P(Xk\Zk)  =  ^pik  xJt>  PjJ.  (117) 

h 

•  After  a  successful  track  extraction  the  sequential  likelihood  ratio  test  is  restarted  and  exploits  the 
remaining  sensor  data  not  used  for  maintaining  the  existent  tracks.  Eventually  other  tracks  can  be 
extracted. 

•  The  sequential  likelihood  ration  test  can  be  used  for  track  censoring:  After  a  decision  in  favor  of  h\ 
we  set:  LR(0)  =  1  and  calculate  LR(k)  from  the  parameters  defining  p(xik\Zli): 

-  Track  confirmation:  LR(k)  >  re-start:  LR(0)  =  1. 

-  Track  Deletion:  LR(k)  <  (eventually  re-extraction  of  the  target) 
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3.6.4  Extension:  Target  Cluster 

The  previously  discussed  track  extraction  scheme  can  directly  be  generalized  to  target  clusters  if  the  number 
n  of  objects  within  the  cluster  is  known.  For  the  sake  of  simplicity  we  let  us  assume  an  ideal  resolution 
capability  of  the  sensor.  In  the  case  of  an  unknown  number  of  objects  we  can  proceed  as  follows: 


1.  Start  the  track  extraction  for  the  target  cluster  with  the  scan  Z\  of  sensor  measurements. 

2.  Assume  a  maximum  number  N  of  objects  within  the  target  cluster,  i.e.  we  have  n  <  N . 

3.  Let  the  a  priori  probability  of  having  n  targets  within  the  cluster  be  given  by:  Pin)  =  jr. 

4.  Consider  hypotheses  hn  assuming  n  individual  objects  within  the  cluster  (ho:  no  object). 

5.  Assume  that  in  the  initial  sensor  data  set  Z\  at  least  one  object  measurement  is  existent. 

1  ^  p(Z*k  |  h  ) 

6.  The  generalized  likelihood  ratio  test  function  is  given  by:  LR(k)  =  —  V  — — - — . 

N  “f  P(Z  \h o) 

7.  The  conditional  likelihoods  p(Zk\h„)  and  p(Zk\ho)  are  iteratively  calculated  as  before. 


3.7  Discussion  of  Examples 

From  our  experiments  with  real  radar  data  we  learned  the  following  lessons  (for  details  see  [21,  18]): 

1 .  IMM-MHT  is  applicable  in  situations  that  are  inaccessible  to  human  radar  operators. 

2.  The  filter  is  rather  robust  and  does  not  critically  depend  on  modeling  parameters  (within  certain  limits). 

3.  Decisive  are  both,  its  multiple  hypothesis  character  allowing  tentative  alternatives  in  critical  situations 
and  the  qualitatively  correct  modeling  of  all  significant  effects. 

4.  Unless  properly  handled,  resolution  conflicts  can  seriously  destabilize  tracking. 

5.  Mono-hypothesis  approximations  to  MHT  (such  as  JPDAF)  are  not  applicable  in  scenarios  as  consid¬ 
ered  in  Figure  3. 

6.  MHT  is  highly  adaptive,  developing  its  multiple  hypothesis  character  only  when  needed. 

7.  Retrodiction  provides  unique  and  accurate  results  from  ambiguous  MHT  output  if  a  small  time  delay 
is  accepted  (some  frames). 

8.  The  maximum  gain  achievable  by  retrodiction  is  roughly  the  same  for  both,  worst-case  modeling  and 
IMM-MHT. 

9.  Algorithms  employing  multiple  dynamics  models  arc  superior  in  that  the  time  delays  involved  arc 
shorter. 

10.  Finally,  it  seems  notable  that  a  very  simplified  modeling  of  the  sensor,  the  target  dynamics,  and  the 
environment  may  provide  reasonable  results  if  applied  to  real  data. 
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